GAMES301课程笔记09-基于调和映射的高质量形变
这个系列是GAMES301-曲面参数化(GAMES 301: Surface Parameterization)的同步课程笔记。课程内容涉及曲面参数化的基本理论与问题描述、面向离散网格的参数化、基于基函数表示的连续参数化、共形参数化、曲面参数化的应用等。本节课主要介绍基于调和映射的网格形变技术。
Deformation
几何形变(deformation)和参数化都可以看做是计算某种映射,它们的主要区别在于参数化一般是高维到低维的映射,而几何形变则主要是考虑同样维度上的映射。

Applications
几何形变在图形学尤其是计算机动画领域有着大量的应用。以传统二维关键帧动画为例,动画的制作过程往往是由原画师先画出角色在关键帧上的形态,这一过程相当于对已有的形状进行变形(shape deformation)。接着动画师会绘制关键帧之间的填充帧,这一过程相当于对形状进行插值(shape interpolation)。几何形变技术可以辅助整个动画制作流程,从而大大减少关键帧动画制作的成本。

因此几何形变技术在工业界有着大量的需求,很多建模动画软件都需要实现相关的功能使用户实现直观、快速且高质量的造型。

Previous Work
目前几何形变技术注意包括基于网格(mesh-based)的方法以及无网格(meshless)方法两大类。基于网格的方法本质是构造分片线性映射,无法保证映射的光滑性;而无网格方法则会直接使用光滑基函数,因此这类方法天生是光滑的。不过现有无网格方法都存在着自身的缺陷,比如说无法保证无翻转单射、无法保证几何扭曲、无法对位置施加约束或者无法保持原形状等。

Planar Map
对于二维平面映射,记\(f: \Omega_1 \rightarrow \Omega_2\)将点\((x, y)\)映射到\((u, v)\)上:
\[\begin{bmatrix} x \\ y \end{bmatrix} \rightarrow \begin{bmatrix} u(x, y) \\ v(x, y) \end{bmatrix}\]对\(f\)的Jacobian矩阵进行SVD可以得到:
\[J = \begin{bmatrix} u_x & u_y \\ v_x & v_y \\ \end{bmatrix} = U \begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \\ \end{bmatrix} V\]其中\(\sigma_1 \geq \sigma_2 \geq 0\)分别表示点\((x, y)\)处两个主方向上的伸缩。
对于二维平面我们还可以使用复数来进行表示,此时\(f\)即为复数域到复数域的映射,其在 \(z=x + i \ y\) 处的导数为:
\[f_z = \frac{1}{2} \bigg( u_x + v_y + i \ (v_x - u_y) \bigg)\] \[f_{\bar{z}} = \frac{1}{2} \bigg( u_x - v_y + i \ (v_x + u_y) \bigg)\]而且这两个导数与奇异值之间存在关系:
\[\sigma_1 = \vert f_z \vert + \vert f_{\bar{z}} \vert\] \[\sigma_2 = \big\vert \vert f_z \vert - \vert f_{\bar{z}} \vert \big\vert\]
利用奇异值或是复数导数可以度量形变的几何扭曲。对于共形映射一般要求两个奇异值尽可能相等,而对于等距映射则会要求两个奇异值都接近于1。

Planar Harmonic Deformation
调和映射是几何形变中非常常用的工具。在二维平面上调和映射满足Laplace方程,它具有很多优良的数学性质包括光滑性(无穷阶可导)、仅在边界上取极值、以及可以通过边界条件唯一确定。

Bounded Distortion Mapping
在设计映射函数时我们一般会要求几何扭曲有界。由于调和映射可以由边界条件唯一确定,因此一个很自然的想法是能否通过控制边界上的几何扭曲来保证全局上的几何扭曲有界。

实际上对于调和映射我们有如下定理保证调和映射在全局上的几何扭曲有界:

同时回忆调和映射可以分解为全纯映射与反全纯映射的和。利用柯西复中心坐标可以对两个映射进行分解,分别用两个cage来控制图形的变化。这样就只需要控制


Bounded Distortion Harmonic Deformation
此时整个计算过程就相当于求解一个约束优化问题。

接下来要求整个映射在全局上扭曲有界,此时需要在边界上进行验证。

基于上面的算法可以得到很高质量的变形效果。



GPU Acceleration
BDHM算法的缺陷在于需要用户指定扭曲的上界,当约束上界过小时可能会出现无解的情况;控制点的变形是通过软约束的形式来施加的,可能无法严格满足;更重要的是它需要求解非凸优化因此很难做到实时。

后续的研究指出可以通过修改优化函数的方式结合GPU加速实现实时形变。这里主要的改进是把ARAP能量改成了对称Dirichlet能量,这样可以把原始的约束优化问题转换成无约束优化问题,从而通过牛顿法进行求解。

和基于网格的分片映射相比,基于调和映射的无网格方法的关键是利用控制点的形变来计算映射函数,相当于使用了边界元。此时的优化问题本质是求解一个小规模的稠密线性方程组,因此适合GPU进行加速。

Newton’s Method
牛顿法是求解非线性无约束优化的基本方法,它通过对优化函数进行二阶泰勒展开来计算最优更新步长。和一阶方法相比牛顿法有二次的收敛速度,可以极大地加速收敛过程。不过需要注意的是二次收敛性依赖于凸的目标函数,当目标函数非凸时牛顿法可能无法收敛。

Isometric Energy
我们希望变形过程接近于刚体变换,因此变形能量需要在Jacobian矩阵的两个奇异值为1时取最小值0。同时变形能量还需要防止图形出现翻转等问题,最好还是连续可微的函数。

目前最常用的变形能量是对称Dirichlet能量。计算几何映射的问题相当于最小化Dirichlet能量在整个区域上的积分,同时满足Jacobian矩阵的行列式大于0。

最后再结合调和映射的性质将全局上的积分转换为边界上的积分,并且使用牛顿法进行求解即可。

Per-Element SPD Hessian
进行求解时的另一个技巧是把Hessian矩阵投影到半正定的空间上,而为了进一步简化计算还可以对控制点再进行抽样来降低计算复杂度。


而在实际使用牛顿法进行计算时可以使用cuBLAS和cuSolver等开源GPU计算库。


Local Injectivity with Lipschitz Continuity
如果再引入Lipschitz连续的条件还可以保证变形在边界上没有翻转。

改进后的BDHM算法基本可以实现实时的几何变形。




总结一下,基于调和映射的二维几何形变有着非常好的数学性质,结合GPU加速的牛顿法可以实现实时的交互变形。

Harmonic Volumetric Deformation
几何形变同样可以应用在三维模型上,此时需要在三维的稀疏点集上计算调和映射。


Volumetric Mapping
对于三维的几何映射,我们需要对二维情况下的一些概念进行推广。

Locally Injective Certification
为了保证映射不出现翻转,我们需要通过采样的方式在采样点上要求无翻转,然后结合Lipschitz连续控制整个区域上的行为。

Challenges in 3D
三维空间中计算几何映射的难点之一在于Jacobian矩阵的最小特征值是没有解析形式的,同时我们也没有相关的定理来将全局的几何扭曲转换为边界上的几何扭曲。

实际上可以构造出边界上几何扭曲有界但内部扭曲大于该上界的调和映射。

Generalization to 3D
因此对于三维空间中的调和映射我们只能通过对内部进行采样的方式来估计和控制几何扭曲。理论上只要采样点足够密集就可以保证整体的几何扭曲有界,当然在这种情况下基本无法实现实时的变形。


Lipschitz Constant for \(\sigma_3\)
对于Jacobian矩阵的最小奇异值\(\sigma_3\),可以使用Hessian矩阵来估计其Lipschitz常数。



实际上当我们使用更高阶的导数时可以得到更紧的上界。

此外还可以把Hessian矩阵的范数加入到目标函数中来得到更光滑的映射。

Optimization
整个优化过程仍然可以使用牛顿法进行计算,也同样可以使用对Hessian矩阵进行投影的技巧。

还需要注意的是几何形变是通过包裹模型的cage顶点来控制的,对于三维的情况控制点往往会远多于二维的情况使得计算上更加复杂。这里引入了降维的技术来减少自由度并提高计算效率。


总结一下,对于三维的体网格同样可以使用调和映射来进行低扭曲高质量的几何变形。
