SIGGRAPH2007的一篇文章翻译:在GPU上直接处理细分曲面(Direct Manipulation of Subdivision Surfaces on GPUs)

在GPU上直接处理细分曲面
(Direct Manipulation of Subdivision Surfaces on GPUs)

Kun Zhou　　 Xin Huang　　 Weiwei Xu 　　 Baining Guo 　 　Heung-Yeung Shum
Microsoft Research Asia


(a)
(b)
图1：细分曲面形变．
(a) 带有位移映射贴图的细分曲面形变，它的控制网格，光滑网格和位移贴图在图2展示．
(b) 带有几何纹理的细分曲面形变．


摘要

我们介绍了一种针对细分曲面形变的算法，适合包括带有位移映射的细分曲面和带有几何纹理的细分曲面．我们的系统让用户自由选择曲面点作为手柄（handle）直接处理曲面．在形变过程中，控制网格的顶点自动调整以使形变曲面满足手柄位置约束．同时保持原曲面的形状和细节．为了最好地保持曲面细节，我们研究出一种梯度域方法能够结合手柄位置约束和细节保持目的到形变能量．对于位移映射细分曲面和带有几何纹理的细分曲面，形变能量是高度非线性和难以处理的在现有的迭代解决方案．为了解决这个问题，我们介绍了一种外壳形变方法，它通过两个稳定的网格形变操作代替每一个不稳定的数字上的迭代步骤．我们的形变算法近使用局部操作，因此对于GPU实现来说是稳定的．结果是比艺术级的多栅格网格形变方案更快的一个实时数量级次序形变系统（The result is a real-time deformation system running orders of magnitude faster than the state-of-the-art multigrid mesh deformation solver）．我们通过一系列例子论证了我们的方法，包括通过操作带有运动捕抓数据的细分曲面实时创作可视的令人满意的角色动画．

关键字：细分曲面，细节保持，位移映射(displacement mapping)，几何纹理

1 引言

细分曲面已经广泛应用到视频创作，商业建模和游戏引擎[DeRose等1998; Warren和Weimer 2002]．通过细分构造曲面给计算机图形工作者带来许多问题，例如任意的拓扑，缩放性，一致性呈现，数值稳定和编码简单[Zorin et al. 2000]．传统的细分曲面主要适合分段光滑建模；[Lee等2000] 介绍了位移映射细分曲面以提高表现效果，该方法整合位移映射[Cook 1984]到细分框架．最近，研究人员更进一步地添加几何纹理到细分表面以使它们成为真正强大的带有复杂细节的表面建模工具．[Peng等2004; Porumbescu 等2005]．

在这篇文章，我们介绍了一种针对细分表面交互形变的算法．我们的算法有下面这些特点：
⋄ 直接处理：用户自由地在目的曲面上选择手柄以代替使用控制网格进行直接处理．为了使曲面形变，用户简单地拖曳手柄到新的位置，我们的算法自动调整控制网格适应手柄位置约束．
⋄ 细节保持：我们的算法能够在产生可视的令人的满意的形变同时有效保持曲面细节．
⋄ 实时表现：我们的算法能够在GPU上实现，对于不同大小的细分曲面能够有100 FPS以上的性能表现．

保持细节对于细分曲面形变是十分重要的．如果没有保持细节，形变曲面会表现出如图10所示的严重扭曲．这促使我们研究一种梯度域的细分曲面形变算法．梯度域方法最近引入到网格形变和编辑[Alexa 2003; Yu等2004; Sorkine等2004]，由于它们能够任意保持曲面细节和产生可视的令人满意的效果是相当著名的．

对于针对细分曲面的梯度域算法一个直接问题是维持细分曲面展现（An immediate issue with a gradient domain algorithm for subdivision surfaces is that of maintaining the subdivision surface representation）．和现有的梯度域方法相似，我们希望直接操作细分曲面并且保持细节，与现有近生成一个形变网格方法不同的是，我们需要一个新的细分控制网格以保证形变结果真正是一个细分曲面．我们通过从曲面网格投影形变能量到控制网格达到这个目的，使用了细分细节函数(subdivision detail function)从控制网格顶点确定曲面网格顶点．

一个更具挑战性的问题是形变过程中保持曲面细节，对于一个没有带有位移映射和几何纹理的细分曲面，细分细节函数是一个简单的线性函数由细分矩阵定义[Warren和Weimer 2002]．在这种情况下形变能量是非线性的并且能够通过高斯－牛顿迭代方法最小化．这和子空间形变技术很相似[Huang等2006]，子空间形变技术使用平均数插值[Ju等2005]以获得稳定快速解决方案．

对于位移细分曲面和带有几何纹理的细分曲面，细分细节函数是非线性的．这导致高度非线性形变能量和高斯－牛顿迭代不在收敛[Huang等2006]．为了处理高度非线性能量，我们介绍了一种外壳形变解决方案．通过使用细分曲面网格作为光滑网格，在它上面叠加位移映射或者几何纹理产生细节网格以创建位移细分曲面和带有几何纹理细分曲面．光滑网格和细节网格分别来源于外壳的内边界和外边界．我们的外壳形变方案在外壳内操作，用两个稳定的形变操作代替数值上不稳定的高斯－牛顿迭代操作：一个针对优化光滑网格，另外一个针对细节网格．通过交替地优化光滑网格和细节网格，我们的方法本质上使用光滑网格形变来计算一个形变能量的高度非线性部分的好的初始值，并且使它容易处理．

我们的算法能够在GPU上实现带来实时表现．关于我们算法的一个关键结论是完全通过局部操作来设计以适应GPU实现．为了平衡CPU和GPU负载，并且利用GPU并行流处理的优点，我们组织细分曲面网格在纹理上如[Shiue et al. 2005]．我们同样预计算需要的逆矩阵，并且把结果加载到GPU作为纹理图象．这样剩下整个迭代在GPU解决带来搞的性能．我们的GPU实现以数量级运行，比使用多栅格的艺术级的形变方案要快．（Our GPU implementation runs orders of magnitude faster than the state-of-the-art fast deformation solver using multigrids [Shi等2006]）

伴随我们介绍的算法，能够在GPU获得高性能高质量形变结果．图1提供了我们算法的形变例子．我们将要证明我们的方法通过更多的例子，我们同样会展示我们GPU形变算法能够使一个动画创作者从静态的细分曲面和动作捕抓数据创作可视的令人满意的实时动画．


2 Related Work
Freeform deformation (FFD) [Sederberg and Parry 1986] embeds
an object inside a volume lattice. The user deforms the object by
manipulating the lattice points. Several extensions have been proposed
to provide a more intuitive user interface by directly manipulating
points [Hsu et al. 1992] or curves [Singh and Fiume 1998]
on the object surface. A recent approach [Botsch and Kobbelt
2005] uses volume-based radial basis functions to deform the object.
Real-time performance on large meshes has been achieved for
deformation with predefined handles.

Figure 2: (a) Control mesh shown in blue. (b) Smooth mesh. (c) Displace-
ment map. See Fig. 1 for the detail mesh of the subdivision surface.


Energy minimization has long been used to deform smooth surfaces
[Welch andWitkin 1992; Botsch and Kobbelt 2004]. [Boier-Martin
et al. 2004] introduces a variational approach to deform subdivision
surfaces. To preserve surface details, they optimize the energy of a
deformation vector field instead of the deformation energy of vertex
positions. With their technique the deformation result is always
a fine mesh with a deformation vector associated with each vertex. This can be quite inconvenient when working with subdivision
surfaces without displacements. Furthermore, since displacements
at vertices are not texture-mapped from a displacement map, this
approach does not scale up well as the subdivision level increases.
Most importantly, the above technique cannot handle geometry textures
that are not displacement maps.
Since the introduction of hierarchical B-Spline editing [Forsey and
Bartels 1988], multiresolution mesh editing techniques [Zorin et al.
1997; Kobbelt et al. 1998; Guskov et al. 2000] have been developed
for detail-preserving deformations by decomposing a mesh
into several frequency bands. A deformed mesh is obtained by first
manipulating the low-frequency mesh and later adding back the
high frequency details as displacement vectors. Recently, [Marinov
et al. 2007] mapped a two-band multiresolution deformation
framework to the GPU. These methods do not support direct manipulation
of the original surface. Also, the displacement vectors
are inserted back independently at each vertex. As a result, artifacts
can appear in highly deformed regions because details are not
coupled and preserved uniformly over the surface.
Gradient domain mesh deformation techniques [Alexa 2003; Yu
et al. 2004; Sorkine et al. 2004; Sheffer and Kraevoy 2004; Zhou
et al. 2005; Lipman et al. 2005; Nealen et al. 2005; Zayer et al.
2005; Au et al. 2006; Huang et al. 2006; Lipman et al. 2006] cast
deformation as an energy minimization problem. The energy function
incorporates position constraints as well as terms for detail
preservation. Minimization of this energy distributes errors globally
over the entire mesh and thus leads to high quality deformation
results. The user can directly manipulate the mesh surface and use
the region of interest to control the scale of manipulation.
Our algorithm combines the strengths of gradient domain techniques
and subdivision surfaces to achieve visually pleasing deformation
and high performance. Recently, [Shi et al. 2006] presents
a fast multigrid solver for gradient domain mesh deformation. Unfortunately,
their GPU implementation does not run much faster
than the CPU version due to the unstructured nature of a general
mesh. Thanks to the regular connectivity and locality-preserving
data access of subdivision surfaces, our deformation algorithm can
be efficiently implemented on the GPU, resulting in a real-time system
which runs orders of magnitude faster than the multigrid solver
of [Shi et al. 2006]. We deem our algorithm a nice complement
to existing GPU-based subdivision techniques [Bolz and Schr ¨oder
2004; Shiue et al. 2005].
Deformation is an active research area and the above review only
summarizes techniques most relevant to our work. Other deformation
approaches include example-based mesh deformation [Sumner
et al. 2005; Der et al. 2006], vector field based shape deformation
[von Funck et al. 2006], and volumetric prism based deformation
[Botsch et al. 2006].


3 细分曲面形变

在这篇文章，一个三角网格 由一个元组 表示，K是一个抽象的单纯复形包括网格的连接信息， 是一个3m维向量，每一个 代表一个顶点位置．


3.1 拉普拉斯形变

我们首先随着[Sorkine等2004;Yu等2004]关于网格的拉普拉斯曲面编辑方案得到一个细分曲面直接处理的公式．细分曲面的控制网格 带有顶点集 ，如图2所示．从控制网格 通过细分 到一定期望层次得到一个带有顶点集 的光滑网格 ．接着在光滑网格 上应用位移贴图[Lee等2000]或细节纹理[Peng等2004]生成带有顶点集 的细节网格 ．

细节网格 是我们希望通过直接处理形变得到的细分曲面．我们通过最小化下面的能量函数进行 的拉普拉斯形变：

　　　　　　　　　　　　　　(1)

其中 是 的拉普拉斯操作矩阵， 是 的拉普拉斯坐标， 是位置约束矩阵， 是约束顶点(直接处理中的顶点)的目标位置． 由于包括局部旋转是顶点位置的非线性函数．我们必须保证形变结果仍然是细分曲面．为此，我们根据控制网格顶点 重写方程(1)．通过细分，位移映射和几何纹理映射，细节网格 的顶点可以从 如下计算：


其中函数 由细分规则，位移映射和纹理映射规则决定．我们称 是细分细节函数(subdivision detail function)．对于一个不带有位移映射和几何纹理的细分曲面， 是由细分矩阵定义的简单的线性函数．但普遍来说 是一个复杂非线性函数．通过用 代替 ，方程(1)可以改写为如下：

　　　　　　　　　　(2)

这是细分曲面拉普拉斯形变的基本公式．首先为控制网格 的顶点集 计算方程(2)，然后应用细分规则，位移映射和几何纹理以达到形变细节网格 ．图3提供一个通过直接处理的细分曲面形变例子．接下来我们介绍如何处理方程2，首先以简单的例子开始，这个例子中位移映射和几何纹理暂不考虑， 是简单的光滑网格 ．

3.2 光滑曲面形变

上一节提到位移映射和几何纹理暂不考虑， 以光滑网格 代替，因此这里简单地把细分细节函数 作为一个线性函数，让 为控制网格 的第 层细分曲面网格，并且 ， 细分矩阵是 ，其中 是 的细分矩阵．细分细节函数 是由细分矩阵定义的线性函数．伴随着 是一个线性函数，那么方程（2）变为：


把系数提出来可以简写为下面的形式：
　　　　　　　　　　　　　　　(3)
，

其中 是光滑曲面 的拉普拉斯操作矩阵，接着我们考察 ．这里我们用 代替 强调它和 的关系．注意到由于 与[Huang等2006]一样是非线性的，所以 是 的非线性函数，方程(3)能够通过不确定高斯－牛顿方法[Steihaug 1995]解决．
　　　　　　　　　 　　　　　(4)

在每次迭代中， 是已知的并且方程(4)能够以最小平方问题(a least squares problem)解决．
　　　　　　　　　　　　(5)

3.3 外壳形变

针对光滑曲面形变的不精确高斯－牛顿方法本质上通过下面线性化：

　　　(6)


其中 是 的雅可比行列式．这个近似值是精确的当 或者步大小 是非常小的．实际上，步大小不总是小的，因为在迭代过程开始时大的步大小经常需要以便快速收敛．幸运的是 确实存在，由于细分细节函数 是线性的并且 的非线性由于非线性拉普拉斯坐标 是单独的．在这种情况下， 仅是适度非线性并且上述的一步线性化方法已经足够．我们通过实验 是在 的范围内．


通常情况下，细节网格 和光滑网格 是不同的并且细分细节函数 是非线性的．这个非线性 导致高度非线性形变能量函数以致不能够使用上述一步线性化方法最小化．为了解决这样的形变能量函数，我们研究出外壳形变方案．为了简单化，我们首先介绍位移细分曲面的外壳形变．假设细节网格 是通过应用位移贴图到光滑网格 创建的，每个顶点 上的顶点 沿着法向量 位移一段距离 ．细节网格的顶点位置 可以如下计算：


是顶点位置 的非线性函数．既然 ， 是 的非线性函数， 是一个 对角矩阵， ．使用 ,我们能够从 计算 如下：



现在细分细节函数 是非线性的由于非线性函数 ．伴随新的 ，方程(2)可以改写为如下：


展开得：


把系数提出来可以简写为下面的形式：
　　　　　　　　　　　　　　　　(7)


是细节网格 的拉普拉斯操作矩阵，我们使用 代替 强调它和 的关系．函数 是高度非线性的由于 和 的非线性．在这样的情况下，[Huang等2006]的一步线性化方法不在足够和适合牛顿－高斯方案通常带来如图4所示的收敛问题．


外壳形变方案是针对方程(7)的迭代方案．从一个带有内边界 和外边界 的薄外壳得到光滑网格 和细节网格 ．在每次迭代中，外壳形变方案优化 和 ．一个迭代从内边界 形变开始，使用位置约束在外边界 上，通过方程(4)根据在 上推断的位置约束实现．推断的约束量由给定的 的上的约束量 导出．光滑网格 形变后， 用于计算方程(7)中 中的非线性拉普拉斯坐标 和位移法向量 ．最后，对于当前迭代细节网格 的形变通过 计算．

特别地，在每次迭代 ，我们首先使用方程(5)计算控制网格 初始估算值并且得到：
　　　　　　　　　　　　　　(8)

其中带有 的 代表光滑网格 上推断的位置约束． 从源位置约束 如下推断．假设一个在细节网格 上的顶点 是根据 受 约束移动的．然后 在光滑网格 上的对应顶点应该根据 受到相同数量 约束移动．

为了计算当前迭代 的形变，我们通过解决下面式子使用已经形变的光滑网格控制顶点 计算 ：

结果如下：
　　　　　　　　　　　　　　　(9)


(a) [Huang等2006]中的一步线性化方案　　　　　　　　　 (b)我们的外壳形变方案
图4：[Huang等2006]中的一步线性化方法和我们外壳形变方案收敛．水平轴代表时间，红色曲线指示形变能量而蓝色曲线指示迭代步的大小，对于动画形变序列参考本文附带演示视频


图3展示了位移细分曲面形变结果．用户直接处理各点在细节网格上．控制网格自动调整并且表面细节精细地保持．图6论证了保持几何细节的重要性．恐龙模型是一个从原始扫描模型通过[Lee等2000]描述算法创建的位移细分曲面．图6(c)是每保持位移贴图细节的形变结果．这是由首先使用3.2节描述算法形变光滑网格然后应用位移贴图得到的的．图6(d)是外壳形变方案的结果．它保持了细节网格的几何细节．就像我们能从图6(e)和(f)中放大版本所看到的．图6(c)中的几何细节是高度压缩的相对于图6(d)．图10提供另外的例子证明了保持细节的重要．

外壳形变方案的稳定性的完整分析已经超过了本文的篇幅．方案后的直觉是不难理解的．外壳形变方案本质上使用光滑网格形变去计算一个 的初始估算，因此让高度非线性的 容易处理．就如[Steihaug 1995]所注明的，高斯－牛顿方法的数值稳定性高度依赖于 非线性． 包括两个非线性部分， 和 ， 更进一步依赖于 ．这复杂非线性让一步线性化方法数值上不稳定，即使用更小的步大小．外壳形变方案用两个数值稳定形变操作取代每一个高斯－牛顿迭代．光滑网格 的形变是稳定的，由于它近涉及 的非线性，因此能够用一步线性化处理．细节网格 的形变是稳定的由于已经形变的 提供了 和 好的初始估算．


我们已经验证了数值稳定性通过一系列实验．图4毕竟了我们外壳形变方法与高斯－牛顿方法的稳定性．就如我们所见，我们的方法更快收敛，而高斯－牛顿方法带有震动地分散．注意到这，我们没有和子空间形变方法[Huang等2006]做一个大体对比．我们仅仅把外壳形变方法和一步线性化方法进行了比较．由于这个原因，我们没有使用平均数插值在这些对比中．代替的是我们调整一步线性化方法根据方程(4)以最小化细分曲面形变能量．既然外壳形变和一步线性化方法使用相同的子空间(例如同样的控制网格和细分方案)，所以这是一个由意义的对比．


图5：带有复杂几何纹理的细分曲面形变(a)光滑网格 (b)几何纹理 (c) 细节网格(d)一个形变结果

我们的算法能够扩展支持位移贴图沿着任意的方向，尽管仅仅沿着法向量方向位移实现了当前．对于普遍的位移，光滑网格的每个顶点的位移方向以一个向量表示在顶点法向量和切向量定义的局部帧．这些局部向量可以存储在额外的纹理．在运行时，我们简单通过局部位移方向和局部帧计算全局位移方向，并且满足外壳形变的结果．


3.4 几何纹理处理

外壳形变方案同样能够处理带有几何纹理的细分曲面．图5展示了一个龙的模型映射了一个鳞状几何纹理．我们使用外壳映射 [Porumbescu等2005]以映射几何纹理到外壳空间覆盖光滑网格 ．首先我们构造一个外壳空间覆盖光滑网格 ．带有偏移的网格 ，它与 具有相同的顶点数量和相同的顶点邻接关系，通过[Cohen等1996]提出的envelope生成算法创建．与位移映射类似， 的每个顶点 vi沿着 处的法向量方向移动一段距离 ．所以偏移网格 的顶点位置能够如下表达：


其中 ，每个 时一个单元向量， 是一个带有 的对角矩阵．

我们定义一个外壳映射通过分解外壳空间( 和 之间的空间)和纹理空间到两个相应四面体的集合．外壳映射由相应四面体的重心坐标定义．给定一个纹理空间的点，我们能够简单定位到它属于的四面体并且计算它的重心坐标．它在外壳空间的相应点定位到带有相同重心坐标的相应的四面体．


伴随着外壳映射，细节网格 的顶点位置能够表达为 和 的线性组合：

其中 是重心坐标矩阵．用 代替 我们得到：



细分细节函数 本质上有着与位移细分表面相同的形式．细节网格 能够通过带有新的矩阵 和非线性函数 的方程(7)形变．
，

再一次， 的非线性结构对于位移映射和几何纹理都是一样的．结果，外壳形变方案同样也能够应用到几何纹理上．注意到我们能够使用其它算法如[Peng等2004]以构造偏移曲面 ．在[Peng等2004]工作中，光滑网格顶点到偏移网格顶点的位移能够任意的．我们能够通过顶点法向量和切向量定义的局部帧来表现．把结果存储为一张纹理．这和我们之前提到的带有任意位移方向的位移映射是相同的．


3.5 细节保持实现

图6：在位移细分曲面保持细节．就如我们能从放大版本里看到的，(c)中的几何细节与(d)相比是高度压缩的
(a)光滑网格
(b)细节网格
(c)没有细节保持的形变结果
(d)带有细节保持的形变结果
(e)图 (c)的放大版本
(f)图 (d)的放大版本


拉普拉斯操作矩阵L能够通过使用[Desbrun等1999]介绍的余切形式构造．对于拉普拉斯坐标 ，我们使用[Huang等2006]介绍的旋转不变表示(rotation-invariant representation)．给定一个没形变网格上的内部顶点 ，它的一圈顶点 和附带的三角形 ，它形变前的拉普拉斯坐标 首先使用 计算．既然拉普拉斯算子是curvature normal的离散近似，它坐落在由附带三角形的法向量跨越(span)的线性空间．接着计算一系列系数 所以：

　　　　　　　　　　　　　　　　
,
要注意的是 和 是为没形变网格预计算的．在每次高斯－牛顿迭代，我们需要计算 ,这里我们用 作为例子解析我们是如何做的．我们首先计算顶点位置 ，然后在每次迭代 使用 计算拉普拉斯值：

　　　　　　　　　　　　(10)

然后我们缩放ei(Vkb )的数量级以保持形变前原来的拉普拉斯值的长度：

　　 　　　　　　　　　　　　　　　(11)

其中 是原始的拉普拉斯值长度．我们当前的系统使用Loop细分方案[Loop 1987]．

4 在GPU上实时形变

我们形变算法的主要部分由各个局部操作组成如细分和拉普拉斯坐标计算以及矩阵向量乘法．这些操作能够在可编程图形硬件上有效实现．接下来我们使用光滑网格形变管线(3.2节)作为例子介绍我们是如何做的．

根据方程(5)，在每个迭代我们需要计算 ，我们在CPU上预计算 和 并且加载结果到GPU上作为两张纹理．可选的，我们可以预计算 并且加载为一张纹理，但是我们没有选择这样做主要由于如下这些原因：首先，因为A是一个稀疏矩阵，计算 仅涉及系数矩阵向量乘法，这并不昂贵．其次，我们希望保持预计算时间足够段以满足用户交互，但计算 会需要更多的预计算时间由于额外的一个稠密矩阵 和一个稀疏矩阵 之间乘法．更少的预计算时间表明能够给予用户更快的交互响应，因为每次用户选择一个新的处理手柄集合我们需要重复预计算过程．有鉴于此，注意到A由两部分组成， 和 ．首部分是通常确定的由于它仅依赖于没形变控制网格和细分层次，第二部分依赖于用户选择的处理手柄，这些手柄在形变中经常改变．所以预计算 提供一个更快的响应给用户．

的计算由两部分组成： 和 ． 是容易计算的因为它直接来源于用户输入，为了计算 ，我们需要通过细分获得光滑网格顶点 ，它能通过[Shiue等2005]介绍的细分核心(subdivision kernel)有效计算．控制网格首先被处理为一系列片段网格(fragment meshes)．共享相同查找表的片段网格被放置到同一组并存储为二维纹理使用螺旋枚举．组中的每个片段网格映射到二维纹理的一行．然后在象素程序(fragment shader)中，为每行用查找表填充必要的细分模版．细分结果(顶点位置和法向量)同样存储到二维纹理用以下一个工序处理或者送入象素缓冲对象(PBOs)用以绘制．一旦我们得到 ，就能根据方程(10)和方程(11)计算 ．

对于最终的 计算，我们执行一次 和 的系数矩阵乘法通过[Bolz等2003]描述的方法，接着进行类似[Kruger和Westermann 2003]所描述的 和 之间的稠密矩阵乘法．