目前已经在UE插件中实现了GPU版本的Marching Cube算法, 接下来继续攻克Dual Contouring算法, 以本文解读又一篇在网格生成领域十分经典的论文《Dual Contouring of Hermite Data》. 相比较Dual Contouring算法, Marching Cube算法的一个比较大的缺点是难以表达尖锐的特征.
参考材料
1. Ju T, Losasso F, Schaefer S, et al. Dual contouring of hermite data[C]//Proceedings of the 29th annual conference on Computer graphics and interactive techniques. 2002: 339-346.
2. Dual Contour通俗解释
3. Dual Contouring Tutorial
4. Isosurface Extraction
5. Givens 旋轉於 QR 分解的應用
6. slides_ch03_ls01
7. 数值稳定性
8. 如何理解施密特（Schmidt）正交化
9. SPACES OF POLYTOPES AND COBORDISM OF QUASITORIC MANIFOLDS
10. cgal-discuss – Re: [cgal-discuss] Difference between Manifold and Quasi-manifold
1. 均匀网格上的Dual Contouring
尽管作者的最终目标是开发一种适用于八叉树的简单等值面方法, 但作者首先探讨了适用于带符号的均匀网格(类似于符号距离场) 的多种等值面技术. 图1左上角展示了带符号的均匀网格的典型示例. 网格中发生符号变化的边会被Hermite数据(亦称为法向量或梯度) 标记, 这些数据包含精确的交点坐标及对应的法向量. 这些Hermite数据可通过轮廓的隐式定义直接计算, 或通过扫描转换封闭多边形网格获得.
1.1 早期等值面方法
基于立方体的方法(如Marching Cubes(MC) 算法及其变体) 会为网格中每个与等值面相交的立方体生成一个或多个多边形. 通常, 这些方法会为等值面与某个特定立方体相交的部分生成一个对应多边形, 其顶点位于等值面与立方体边的交点处. 图1右上角展示了由左侧带符号的网格生成的MC等值面的二维示例. 图2左侧显示了一个三维示例: 将函数$f(x, y, $$ z) = 1 – x^2 – y^2 – z^2$的零等值面生成为一个球体. 该等值面由一系列多边形组成, 这些多边形近似表示等值面在网格中每个立方体内的局部形态.
对偶方法(如Gibson 1998年提出的SurfaceNets算法) 会为每个与等值面相交的立方体生成一个位于等值面或其附近的顶点. 对于网格中发生符号变化的每条边, 与包含该边的四个立方体相关联的顶点会被连接起来, 形成一个四边形. 最终生成一个连续的多边形表面以逼近等值面. 图2右侧展示了使用SurfaceNets方法生成的同一球体示例. 需要注意的是, SurfaceNets方法生成的多边形网格在标准拓扑意义上与MC方法生成的网格对偶: SurfaceNets网格的顶点对应MC网格的面, 反之亦然. 对偶方法通常能生成具有更优纵横比的多边形网格, 因为网格顶点可以自由地在立方体内移动, 而不像基于立方体的方法那样受限于网格边.
$\\$ 扩展的Marching Cubes(Extended Marching Cubes, EMC) 方法是一种结合基于立方体方法和对偶方法的混合算法. EMC通过分析立方体边上的交点法向量来检测立方体内是否存在尖锐特征. 若某个立方体的法向量位于用户指定的锥形区域内, 则判定该立方体无尖锐特征. 此时, EMC使用标准Marching Cubes(MC) 方法生成一个或多个多边形. 对于包含尖锐特征的立方体, 该方法会生成一个顶点, 其位置位于一个二次函数的极小化点处, 该二次函数定义如下所示:$\begin{align}E[x] = \sum_{i} (n_i \cdot (x – p_i))^2, \label{QEF} \end{align}$其中, $(p_i, n_i)$对应等值面与立方体的交点(及其单位法向量). 一旦该顶点定位完成, 该方法会通过生成三角扇(指以单一顶点为中心辐射连接多边形顶点的三角化方式) 连接到立方体边界上的边. 最终, 若两个相邻立方体均包含尖锐特征顶点, 则三角扇在它们的公共面上生成的三角形对会翻转其公共边, 形成尖锐特征边. 图1左下角展示了由EMC生成的等值面二维示例.
1.2 基于Hermite数据的Dual Contouring算法
EMC方法的主要优势在于其利用Hermite数据和二次误差函数(QEF) 来定位包含几何特征的立方体顶点. 这种Hermite方法能够生成同时包含尖锐顶点和尖锐边的等值面. 然而, 这一方法的缺点是需要显式检测此类特征, 并在这些情况下执行特定的后处理步骤. 作为EMC方法的替代方案, 作者提出以下针对Hermite数据的对偶等值面方法:
$\\$ 1) 对每个发生符号变化的立方体, 生成一个位于方程$\eqref{QEF}$二次函数极小化点处的顶点;
$\\$ 2) 对每条发生符号变化的边, 生成一个四边形, 连接包含该边的四个立方体的极小化顶点.
$\\$ 该方法是EMC方法与SurfaceNets方法的一种有趣混合策略. 它采用EMC方法的特征顶点规则来定位等值面的所有顶点, 同时利用SurfaceNets方法确定这些顶点的连通性(需注意, SurfaceNets方法使用了完全不同的规则来定位等值面上的顶点). 通过使用QEF定位等值面的所有顶点, 该方法避免了显式检测特征的需求. 等值面上的顶点仅需根据Hermite数据对应的法向量一致性进行定位. 图1右下角展示了由左上角Hermite数据生成的对偶等值面的二维示例.
$\\$ 图3展示了通过在$64^3$网格上Dual Contouring算法处理Hermite数据建模的机械零件3D示例. 左图显示了零件的平滑着色版本, 右图展示了Dual Contouring算法生成的多边形网格. 该模型的交点和法向量由封闭细分曲面生成, 标定模型内部 / 外部的符号距离场通过标准扫描转换算法计算(如[Foley et al. 1995] 所述).
1.3 QEF的表示与最小化
在此, 需要对QEF的表示与最小化方法作几点说明. 方程$\eqref{QEF}$中的函数$E[ $$ x]$由一组交点$p_i$和法向量$n_i$构建而成. 该函数$E[x]$可表示为内积形式$( $$ Ax – b)^T (Ax – $$ b)$, 其中矩阵$A$的行向量为法向量$n_i$, 向量$b$的元素为$n_i \cdot $$ p_i$(即法向量与交点的内积). 通常, 二次函数$E[x]$会被展开为以下标准形式:$\begin{align}E[x] = x^T A^T Ax – 2x^T A^T b + b^T b, \end{align}$其中, 矩阵$A^T A$是一个对称的$3 \times 3$矩阵, $A^T b$是一个长度为3的列向量, $b^T b$是一个标量. 这种展开的优势在于, 仅需存储$A^T A$, $A^T b$和$b^T b$(共10个浮点数), 而无需存储完整的矩阵$A$和向量$b$. 此外, 通过求解正规方程$A^T A\hat{x} = A^T b$, 即可计算出$E[x]$的极小值点$\hat{x}$.
$\\$ 这种表示的一个缺点在于其数值不稳定性. 例如, 当从平坦区域采样交点和法向量构建$E[x]$时, 在浮点数运算中计算$E[x]$的值可能会产生显著误差. 其数值不稳定性的来源为:
$\\$ $\cdot$ 法向量共线性: 平坦区域的法向量几乎平行, 导致矩阵$A^T A$接近奇异(行列式趋近于0), 求解正规方程时误差放大.
$\\$ $\cdot$ 大数值(如$b^T b \approx 10^6$) 与小数值(如$E[x] \approx 0$) 的混合运算中, 有效数字被截断, 导致计算结果不可靠.
$\\$ 以$256^3$网格(如图4所示) 为例, $b^T b$的量级可能达到$10^6$级别. 由于浮点数仅能保证六位十进制有效数字的精度, 若在原始平坦区域(理论上$E[x]$应为零) 的点上计算$E[x]$, 其结果误差可能达到1($\approx 10^6 \times 10^{-6}$) 的量级.
解决这一问题的一种可行方案是, 在表示$A^T A$, $A^T b$和$b^T b$时改用双精度浮点数代替单精度浮点数. 通过使用双精度数, 平坦区域示例中的$E[ $$ x]$误差将降低至$10^{-6}$($\approx 10^6 \times 10^{-12}$) 量级. 然而, 采用双精度的主要缺点是存储一个QEF所需的空间将翻倍. 在本文的应用中, 这一方案并不可行, 因为程序性能受存储空间限制(而非时间限制).(详见最后一节的详细说明.)
$\\$ 一种替代的QEF表示方法可在仅使用单精度浮点数的同时达到双精度的准确性, 其基于QR分解[Golub与Van Loan 1989]. 若将矩阵$A$附加列向量$b$后形成的矩阵记为$(A \ | \ b)$, 该分解的核心思想是计算一个正交矩阵$Q$, 使得$Q$与$(A \ | \ b)$的乘积为如下形式的上三角矩阵:$\begin{align} \begin{pmatrix}
x & x & x & x \\
0 & x & x & x \\
0 & 0 & x & x \\
0 & 0 & 0 & x \\
0 & 0 & 0 & 0 \\
\cdots & \cdots & \cdots & \cdots
\end{pmatrix} = \begin{pmatrix}
\hat{A} & \hat{b} \\
0 & r \\
0 & 0 \\
\cdots & \cdots
\end{pmatrix}. \label{upper_triangular_matrix_form_of_augmented_matrix_of_matrix_A} \end{align}$其中, $\hat{A}$是一个上三角$3 \times 3$矩阵, $\hat{b}$是一个长度为3的列向量, $r$是一个标量. 该矩阵$Q$可以表示为一系列Givens旋转的乘积, 其中每个旋转将$(A \ $$ | \ b)$的下三角部分中的单个元素归零.
$\\$ 由于任何正交矩阵$Q$均满足关系式$Q^T Q = I$, 函数$E[x]$可被重写为:$$\begin{aligned} (Ax – b)^T (Ax – b) &= (Ax – b)^T Q^T Q(Ax – b) \\ &= (QAx – Qb)^T (QAx – Qb) \\ &= (\hat{A}x – \hat{b})^T (\hat{A}x – \hat{b}) + r^2. \end{aligned}$$以这种形式计算$E[x]$时, 首先计算向量$\hat{A}x – \hat{b}$与自身的点积, 然后再加$r^2$. 如此一来, 问题转化为关于上三角系统$Rx = c$的求解, 而不涉及病态矩阵$A^T A$; 此外, 一般来说, 上三角矩阵$R$的条件数亦低于矩阵$A$, 降低了数值误差的影响. 回到之前的平坦区域示例, 作者注意到$b$的元素量级在$10^3$左右, 因此当$x$选自平坦区域时, $\hat{A}x – \hat{b}$的元素量级为$10^{-3}$. 这是因为立方体的尺寸限制了残差的理论最大值:
$\\$ $\cdot$ 假设立方体边长为$a$(例如$a = 1$, 归一化单位), 则$|x – p_i| \le s \sqrt{3} \approx $$ 1.732$(单位立方体对角线长度).
$\\$ $\cdot$ 由于$n_i$是单位法向量, 残差的绝对值满足:$$ |\hat{A}x – \hat{b}| = |Q| \sum_i |n_i \cdot (x – p_i)| \le \sum_i |x – p_i| \le \sum_i (s \sqrt{3}) \approx \sum_i 1.732. $$由此, 计算得到的$E[x]$值将处于$10^{-6}$量级.
$\\$ 如果$\hat{A}$是非奇异矩阵, 则可以通过回代法求解$\hat{A}x = \hat{b}$来计算极小值点$\hat{x}$. 然而, 在Dual Contouring算法中, $\hat{A}$通常由近似共面的噪声法向量计算得到. 此时, 矩阵$\hat{A}$接近奇异. 结果, 极小值点$\hat{x}$可能位于定义立方体之外. 为解决这一问题, 作者对$\hat{A}$进行SVD, 并通过截断其较小的奇异值形成伪逆矩阵(方法参考Kobbelt et al. 2001; Lindstrom 2000]). 基于实验验证, 作者通常截断所有奇异值量级小于0.1的部分. 利用所得的伪逆矩阵, 作者近似求解$\hat{A}x = \hat{b}$, 同时最小化$\hat{x}$与所有交点$p_i$构成的多边形的质心之间的距离.
1.4 纹理轮廓线建模
在图2中, 两种等值面算法生成的表面均分割了从负空间(空空间) 到正空间(实体空间) 的过渡. 在真实环境中, 实体并非由单一均匀材质构成. 实际上, 实体由多种材质共同组成, 每种材质会在表面上形成独特的纹理. 图5展示了一个由两种材质组成的立方体示例: 一种是金色材质, 通过将汉字沿立方体方向拉伸贯穿而成; 另一种是红色材质, 构成了立方体的其余部分. 需注意, 金色材质是真正的实体(而非表面纹理), 因为汉字贯穿了整个立方体(如右图中球形切割后的立方体所示).
通过将实体分割为不同材质, 这种划分可以隐式地建模: 用材质索引替代原本表示空空间(负值) 和实体空间(正值) 的符号-和+. 在此表示方法中, 每个网格点均有一个对应特定材质的索引. (参见[Bloomenthal and Ferguson 1995; Bonnell et al. 2000] 了解类似方法的示例.) 图6展示了一个包含三种不同索引的二维网格; 灰色和黑色网格点表示不同的实体材质, 而白色网格点表示空空间. 与之前相同, 存在索引变化的边缘由精确的交点和法向量标记. 在等值面生成过程中, 这些Hermite数据通过方程$\eqref{QEF}$定义与网格中每个立方体相关联的QEF. 随后, 对于每条存在索引变化的边缘, Dual Contouring算法生成一个四边形, 连接包含该边缘的四个立方体的QEF极小值点. 图6左侧展示了分隔三种材质的对偶等值面.
若观察者仅能从空空间视角观察, 可以针对由多种材质组成的实体优化此等值面方法. 特别是, 由固 / 固边缘生成的四边形面片在空空间中不可见(当观察者处于空空间时, 实体内部材料之间的边界(如红 / 金分界) 无法被直接观测). 剩余的四边形面片对应固 / 空边缘, 并可通过边缘实体端点的材质属性进行纹理映射. 在图5中, 底层网格包含三种材质: 空空间, 金色材质和红色材质. 生成的对偶等值面由红色 / 空边缘生成的红色四边形和金色 / 空边缘生成的金色四边形组成. 红色与金色材质之间的边缘(红 / 金边缘) 不生成四边形面片(通过忽略内部边界, 减少网格面数).
$\\$ 当三种或更多材质在单个立方体内交汇时, Dual Contouring算法会将极小值顶点放置在交汇点或其附近. 这种定位方式能够高度还原表面凸起文字和字符的轮廓(图7右侧局部放大图展示了这一效果的示例). 基于立方体的轮廓方法将顶点约束在3D网格的边上, 因此难以实现此类效果.
在多材质场景下, 可以对双材质情况中使用的CSG操作进行多种有趣的扩展. 取代传统的CSG操作, 作者采用单一操作Add, 用新材质覆盖现有模型的一部分. 例如, 图5右上角展示的球形切割等减法操作, 可通过向模型中添加一个空空间球体来实现. Add的一种有用变体是Replace操作, 它仅覆盖模型的实体部分. 该操作还可用于模拟对轮廓部分区域的纹理化. 例如, 图5下方的图像展示了将实体部分替换为蓝色材质后, 再通过球体进行切割的效果.
2. 自适应Dual Contouring算法
之前的Dual Contouring算法存在一个显而易见的缺陷: 其设计基于均匀网格. 在实际应用中, 均匀网格的大部分存储空间用于保存同质立方体(即所有顶点符号相同的立方体), 仅有少量异质立方体(即包含等值面交点的立方体) 参与表面生成. 为避免这种空间浪费, 一种解决方案是用八叉树替代均匀网格. 本节将描述一种基于自适应简化八叉树的Dual Contouring算法, 其中八叉树的叶节点包含QEF. 该方法本质上是[Lindstrom 2000] 提出的均匀简化方法的自适应变体. 作者的方法包含以下三个步骤:
$\\$ 1) 生成带符号的八叉树, 其同质叶节点已最大程度折叠(合并).
$\\$ 2) 构造QEF并简化八叉树: 为每个异质叶节点构建QEF, 并基于这些QEF简化八叉树.
$\\$ 3) 递归生成多边形: 对简化后的八叉树递归生成表面多边形.
$\\$ 需要注意的是, 作者的方法与Lindstrom的方法也有区别: 作者从带符号的八叉树中直接生成多边形, 而非对现有网格中的多边形进行折叠操作.
$\\$ 在第一步中生成带符号的八叉树相对简单. 对于隐式模型或面片模型, 可以通过对模型进行空间划分, 以自上而下的方式递归构建该八叉树. 对于均匀网格上的带符号的数据, 则可以通过自下而上的方式递归合并均匀区域来生成八叉树. 接下来的两个小节将更详细地探讨该流程的第二步和第三步.(需注意, 此处描述的自适应方法无需修改即可适用于多材质场景.)
2.1 基于QEF的八叉树简化
作者对自适应方法第二步的处理是: 使用方程$\eqref{QEF}$为每个异质叶节点构建对应的QEF. 需注意, 该QEF的极小值点对应的残差可估计顶点对原始几何形状的逼近精度[Garland and Heckbert 1998]. 对于八叉树的简化, 作者的方法是: 通过累加当前节点子树中所有叶节点的QEF, 为八叉树内部节点生成新的QEF. 若某内部节点的QEF残差小于给定容差, 则将其折叠为叶节点.
$\\$ 作者对这一方法所做的唯一修改是修改QEF的内部表示形式, 以利用前文讨论的QR分解技术. 为此, 作者通过方程$\eqref{upper_triangular_matrix_form_of_augmented_matrix_of_matrix_A}$中的上三角矩阵元素, 用10个浮点数表示一个QEF. 两个QEF的合并操作对应于将其各自的两个上三角矩阵进行合并, 形成一个$8 \times 4$的联合矩阵, 其形式如下:$$\begin{pmatrix}
x & x & x & x \\
0 & x & x & x \\
0 & 0 & x & x \\
0 & 0 & 0 & x \\
x & x & x & x \\
0 & x & x & x \\
0 & 0 & x & x \\
0 & 0 & 0 & x
\end{pmatrix}.$$随后, 通过一系列Givens旋转将矩阵重新转换为方程$\eqref{upper_triangular_matrix_form_of_augmented_matrix_of_matrix_A}$的上三角形式. 由于Givens旋转的正交性, 合并系统的QEF等于未合并系统QEF的和. 需注意, 将合并后的矩阵重新转换为上三角形式的速度比标准表示法中直接相加10个浮点数更慢(约需150次算术运算). 然而, 这种表示法的数值稳定性提升带来了更优的减面效果.
$\\$ 图8直观展示了QR表示法稳定性的优势. 该图中的网格展示了对图4中寺庙模型在误差容限为0.014时的两种减面结果: 左侧网格使用QEF的标准表示法计算, 右侧网格则使用QEF的QR表示法计算.(为直观理解规模, 寺庙模型定义在一个$256^3$单元的网格上.) 由于标准表示法的不稳定性引入的数值误差(因残差估计失真, 如误判应合并区域为需保留, 需额外多边形保持几何失真), 左侧网格包含78K个多边形, 而右侧网格仅需36K个多边形.
2.2 关于简化八叉树的多边形生成
基于这一简化后的八叉树, 作者的下一个任务是适当调整Dual Contouring算法的多边形生成阶段. 对于基于立方体的方法, 从八叉树生成等值面的问题已被广泛研究[Bloomenthal 1988; Wilhelms和Gelder 1992; Livnat等1996; Shekhar等1996; Westermann等1999; Frisken等2000; Cignoni等2000]. 通常, 这些方法会限制八叉树的相邻叶节点层级差异最多为1(即”受限” 八叉树), 并通常执行某种类型的裂缝修复以确保等值面闭合. [Perry和Frisken 2001] 描述了一种基于枚举八叉树叶节点相关边的带符号的八叉树SurfaceNets算法变体. 具体而言:
$\\$ $\cdot$ 对于每条表现出符号变化的边, 生成所有连接该边所包含的任意三个不同立方体对应顶点的三角形.
$\\$ 在此过程中, 可能生成冗余三角形, 其来源主要为:
$\\$ $\cdot$ 在八叉树简化中, 叶立方体可能包含多个相邻材质边界(如汉字笔画交叉处).
$\\$ $\cdot$ 直接生成所有可能三角形会导致重叠或冗余(如图5中球形切割的空空间区域误生成三角形).
$\\$ 为避免生成冗余三角形, 该方法根据三角形对应边在叶立方体内的相对位置剔除部分生成的三角形. 在最终网格中, 每条边对应于轮廓上的一个或两个三角形(即一个四边形面片). 该方法的主要缺点(作者已承认) 是偶尔会产生带有裂缝的轮廓. 裂缝的成因主要为:
$\\$ $\cdot$ 受限八叉树的层级差限制(相邻叶节点最多差1级) 可能无法完全覆盖所有拓扑情况.
$\\$ $\cdot$ 剔除冗余三角形后, 部分边可能仅对应单个三角形, 导致网格不连续(如图8左侧标准表示法的78K面因数值误差产生裂缝).
$\\$ 作者通过识别这些导致问题的配置, 并在八叉树上执行额外细分来规避这些问题. 其算法流程大致如下所示:
$\\$ 1) 步骤1: 边相对位置检测.
$\\$ $\quad$ $\cdot$ 对每条边, 计算其在叶立方体内的位置(如靠近立方体中心或顶点).
$\\$ $\quad$ $\cdot$ 若多个三角形共享同一组边且位置重叠, 则剔除冗余三角形(如保留最靠近材质交界的三角形).
$\\$ 2) 步骤2: 四边形面片生成. 对每条边, 生成1或2个三角形.
$\\$ $\quad$ $\cdot$ 若边仅被一个立方体包含, 则生成一个三角形(如平坦表面).
$\\$ $\quad$ $\cdot$ 若边被两个立方体包含(如材质交界), 则生成两个三角形组成四边形面片.
$\\$ 3) 步骤3: 裂缝修复.
$\\$ $\quad$ $\cdot$ 识别导致裂缝的配置(如边仅对应单个三角形且相邻区域未细分).
$\\$ $\quad$ $\cdot$ 对相关叶立方体执行额外细分(增加八叉树层级), 强制生成匹配的三角形对.
$\\$ 作者提出了一种更简单的Dual Contouring算法规则, 用于带符号的八叉树, 从而避免额外细分的需求. 该规则基于以下观察: 只有那些不包含相邻叶立方体边的边才会生成多边形. 作者将这类边称为八叉树的最小边. 因此, 本文的多边形生成规则为:
$\\$ $\cdot$ 对于每条表现出符号变化的最小边, 生成连接包含该边的立方体的极小值顶点的多边形.
$\\$ 这一规则的特性在于, 它能够为任何简化后的八叉树生成闭合的多边形网格. 具体而言, 网格中的每条边均被偶数个多边形包含. 为证明这一特性, 作者观察到对偶等值面上的边由面相邻的叶立方体对生成. 覆盖其共同正方形面边界的最小边始终表现出偶数次符号变化, 因为正方形的边界是一条闭合曲线. 规则总会生成偶数个多边形包含该边. 例如, 在均匀情况下, 共同正方形面由四条连续的边组成. 这四条边链可能表现出两次或四次符号变化, 从而生成两个或四个包含共同边的多边形.(也可构造带符号的八叉树生成包含6个或更多多边形共享同一条边的对偶等值面.) 图9展示了机械零件的三种减面近似结果. 需注意, 最右侧网格经历了拓扑变化.
在八叉树的过渡区域(即一个粗糙网格立方体与四个精细网格立方体面相邻的区域), 该规则会生成三角形而非四边形. 位于公共粗糙网格面中央的最小边仅被三个立方体包含, 这些边会生成三角形, 从而在粗糙网格四边形和精细网格四边形之间形成过渡. 图8展示了通过等值面生成简化八叉树时产生的许多此类过渡三角形的示例.
$\\$ 需注意, Perry / Frisken规则通过枚举八叉树中的边并定位包含该边的立方体来实现. 这一邻接查找过程需要遍历八叉树的上下层级或显式维护相邻立方体之间的链接. 相较于枚举边并主动查找邻接关系, 作者提出了一种递归枚举方法, 用于定位包含共同最小边的所有立方体集合. 为简化说明, 作者以四叉树为例解释该枚举方法, 同时说明其原理可自然扩展至八叉树.
$\\$ 该枚举过程的关键在于两个递归函数faceProc[$q$]和edgeProc[$q_1$,$q_2$].
$\\$ $\cdot$ 给定四叉树中的一个内部节点$q$, faceProc[$q$]会递归调用自身与$q$的四个子节点, 并对$q$的所有四组边相邻子节点对调用edgeProc.
$\\$ $\cdot$ 给定一对边相邻的内部节点$q_1$和$q_2$, edgeProc[$q_1$,$q_2$]会递归调用自身与跨越$q_1$和$q_2$公共边的两组边相邻子节点对.
$\\$ $\cdot$ 图10展示了这两个函数的相互递归结构.
当$q_1$和$q_2$均为四叉树的叶节点时, 对edgeProc[$q_1$,$q_2$]的递归调用终止. 此时, edgeProc[$q_1$,$q_2$]已掌握生成$q_1$和$q_2$共享的最小边相关线段所需的全部信息. 需注意, 该方法的运行时间与四叉树的大小呈线性关系, 因为每个四叉树方块对应一次faceProc调用, 每条边对应一次edgeProc调用.
$\\$ 等值面生成八叉树需要三个函数: cellProc[$q$](单元处理函数), faceProc[$q_1$, $q_2$](面处理函数) 和edgeProc[$q_1$,$q_2$, $q_3$, $q_4$](边处理函数).
$\\$ $\cdot$ cellProc函数会生成8次cellProc调用(递归处理8个单元), 12次faceProc调用(处理12对边相邻子单元), 6次edgeProc调用(处理6对共享公共边的子单元对).
$\\$ $\cdot$ faceProc会生成4次faceProc调用(递归处理4对边相邻子单元), 4次edgeProc调用(处理4对共享公共边的子单元组).
$\\$ $\cdot$ edgeProc会生成2次edgeProc调用(递归处理跨越公共边的两组子单元).
$\\$ $\cdot$ 当所有$q_i$均为八叉树叶节点时, 对edgeProc[$q_1$,$q_2$, $q_3$, $q_4$]的递归调用终止于八叉树的最小边处.
3. 保持拓扑的减面
减面方法如[Rossignac和Borrell 1993; Lindstrom 2000] 的特性是, 多边形网格的拓扑连通性可能在减面过程中发生变化. 更复杂的方法如[Stander和Hart 1997; Gerstner和Pajarola 2000; Wood等2000; Guskov和Wood 2001] 被开发用于在减面过程中保持网格的连通性. 遗憾的是, 在本文的场景中, 不仅轮廓的拓扑连通性, 其纹理区域的连通性也可能变化. 图7的左侧展示了汉字笔画在减面中合并的示例, 导致难以辨认. 尽管这些拓扑变化并非总是不可取, 作者希望在减面过程中能够选择性地保留轮廓及其纹理区域的拓扑连通性. 具体而言, 给定八叉树中一个内部节点(其八个子节点均为叶节点), 需要一种基于这些叶节点角点处的符号(或索引) 的测试方法, 以确保在节点折叠时对偶等值面及其纹理区域的拓扑连通性得以保留.
3.1 双符号情况
考虑一个由八个叶立方体组成的粗糙立方体. 这八个叶立方体的角点符号定义了一个$3 \times 3 \times 3$的网格, 而该网格的角点又构成了一个$2 \times 2 \times 2$的粗糙网格, 如下图所示. 作者的目标是开发一种测试方法, 以判断由该精细网格生成的对偶等值面是否与由粗糙网格生成的对偶等值面拓扑等价.
在提出测试方法之前, 首先回顾: 若一个$d$维等值面在局部与$d$维圆盘拓扑等价, 则其为流形. 由于立方体有12条边, Dual Contouring算法最多可在立方体中心顶点处生成12个多边形. 对于立方体上大多数常见的符号配置, 这些多边形在该顶点处定义了一个流形. 然而, 存在某些符号配置会使对偶等值面成为非流形.(这些配置对应于传统基于立方体的方法中的”歧义” 符号配置.) 基于此定义, 安全性测试包含以下三个检查步骤:
$\\$ 1) 检查粗糙立方体的对偶等值面是否为流形. 若否, 则终止减面.
$\\$ 2) 检查每个精细立方体的对偶等值面是否为流形. 若否, 则终止减面.
$\\$ 3) 检查精细等值面在粗糙立方体的每个子面上是否与粗糙等值面拓扑等价. 若否, 则终止减面; 否则允许安全折叠.
$\\$ 前两个检查步骤将减面过程限制在生成流形对偶等值面的范围内.(需要注意的是, 如果精细叶立方体本身是前一次折叠的结果, 则可以省略第二个检查.) 在实际应用中, 这种限制是可接受的, 因为大多数高分辨率等值面本质上是流形的, 而非流形结构通常由不安全的减面操作导致. 对于前两个检查步骤, [Gerstner和Pajarola 2000] 描述了一种简单测试方法, 用于判断单个立方体的等值面是否为流形. 其核心思想是: 将立方体角点符号相同的边反复折叠为单个顶点. 此时, 该立方体的等值面是流形当且仅当折叠后的结果为单条边. 该测试结果可对所有可能的单个立方体符号配置预先计算, 并存储在一个大小为$2^8$的查找表中.
$\\$ 第三种检查通过以下步骤验证粗糙网格与精细网格等值面的拓扑等价性:
$\\$ 1) 首先检查粗糙立方体边上的拓扑等价性.
$\\$ 2) 接着检查粗糙立方体面上的拓扑等价性.
$\\$ 3) 最后检查粗糙立方体内部的拓扑等价性.
$\\$ 这些检查可通过一系列对$3 \times 3 \times 3$带符号的网格的符号比较实现:
$\\$ $\cdot$ 粗糙边中点的符号必须与该边两个端点中的至少一个符号一致.
$\\$ $\cdot$ 粗糙面中心的符号必须与该面四个角点中的至少一个符号一致.
$\\$ $\cdot$ 粗糙立方体中心的符号必须与该立方体八个角点中的至少一个符号一致.
$\\$ 图11展示了三个带符号的四叉树的减面候选案例. 左侧四叉树的对偶等值面包含两个不同的连通组件. 在此情况下, 第一个检查判定该折叠为不安全操作, 因为折叠后的四叉树等值面为非流形结构. 中间四叉树的对偶等值面同样包含两个连通组件, 但第三个检查判定其不可安全折叠, 因为该四叉树的左侧边无法通过拓扑等价性验证. 右侧四叉树的符号配置满足所有三项检查条件, 因此该四叉树可被安全折叠.
这些符号检查的正确性证明基于按维度递增顺序建立粗糙立方体子面的拓扑等价性. 图9中的右侧网格展示了机械零件减面版本的一个示例, 该减面经历了导致网格断开的拓扑变化. 图9中的中间网格展示了在拓扑检查保护下完成安全减面后的结果, 该检查阻止了进一步的不安全减面.
3.2 多材质情况
前文讨论的等值面生成与减面方法的一个显著优势在于, 这些方法能够轻松处理多材质情况(无需额外复杂度). 幸运的是, 前一小节描述的安全性测试也可通过一个微小调整扩展至多材质区域的等值面.
$\\$ 表面上, 挑战在于多材质区域的等值面在双材质意义下本质上是非流形的. 例如, 图12展示了三个分离三种材料的对偶等值面示例. 其中两个等值面在某个顶点处存在三种材质的交界. 需要注意的是, 若单独考察每种材质区域的边界, 则仍可判定该部分对偶等值面是否为流形. 具体而言, 若每种材质区域的边界均为流形, 则该多材质对偶等值面为Quasi流形(根据Quasi流形的定义, 一对$d$维Cells只能通过共享的$d – 1$维Cell相连, 详见参考材料9, 10). 在双材质情况下, Quasi流形等于流形.
多材质安全性测试判断经简化后的对偶等值面(通常为Quasi流形) 是否是保持拓扑的. 与之前类似, 这一限制并非特别成问题, 因为多材质等值面的大部分区域本质上是Quasi流形. 这一新测试同样分为三个阶段, 与双材质测试完全相同, 但需替换第一阶段和第二阶段的检查规则: 将判断立方体内部等值面是否为流形的测试替换为判断其是否为Quasi流形的等效测试, 而第三阶段的索引测试保持不变.
$\\$ 与流形情况类似, 多材质立方体的Quasi流形测试涉及将立方体中端点具有相同材质索引的边进行折叠. 此时, 该立方体的对偶等值面是Quasi流形当且仅当折叠后的边图形成一个单形(即点, 线段, 三角形或四面体). 与双材质情况类似, 该测试的值可以预先计算并存储在一个大小为$4^8$的查找表中.(若立方体包含5种或更多不同的材质索引, 其边图无法折叠为单形.) 该测试的正确性可通过以下方法验证: 选择立方体上的一个材质索引, 并将所有其它索引视为等效材质. 此时折叠后的边图形成线段, 表明对应所选索引的等值面部分为流形.
$\\$ 图12展示了两个多材质四叉树的简化候选案例. 左侧四叉树因第三阶段检查在底部边失败(如中点符号与端点不一致) 而被拒绝简化. 中间四叉树通过了所有三个检查, 并折叠成右侧的四叉树. 需要注意的是, 简化后的四叉树等值面为Quasi流形, 因为该正方形的边图折叠后形成了一个三角形.
$\\$ 图7展示了图5中汉字立方体的两个减面版本. 左侧版本在减面过程中未采用任何保持拓扑的措施, 导致原本分离的汉字笔画部分融合在一起. 右侧版本通过多材质测试实现了保持拓扑的减面, 保留了字符各区域的独立性(减面后仍保持分离).
附注:
$\\$ 名称中的”对偶” 源于网格中的单元在生成的网格中转换为顶点, 这一特性与对偶图相关.