GAMES301课程笔记04-无翻转参数化方法(初始无翻转)
这个系列是GAMES301-曲面参数化(GAMES 301: Surface Parameterization)的同步课程笔记。课程内容涉及曲面参数化的基本理论与问题描述、面向离散网格的参数化、基于基函数表示的连续参数化、共形参数化、曲面参数化的应用等。本节课主要介绍初始无翻转情况下的的无翻转参数化方法。
计算无翻转参数化的另一种思路是以一个初始无翻转的参数化作为起点，然后通过不断优化的方式来修正几何扭曲。
这类方法一般会以Tutte’s embedding作为初始化，然后在无翻转映射的空间中通过迭代的方式来求解几何扭曲的最小化问题。整个求解流程可以抽象为下图：
为了保证优化过程中不会出现翻转的现象，在整个目标函数中还会加入一个障碍函数(barrier function)。当三角形趋向于退化或是翻转时，障碍函数会快速达到无穷大从而避免这种情况。除了单独设置障碍函数外还可以使用几何扭曲的度量函数来实现类似的功能，比如说MIPS函数就可以避免退化的情况：
\[\frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1} = \frac{\sigma_1^2 + \sigma_2^2}{\sigma_1 \sigma_2} = \frac{\Vert \mathbf{J}_t \Vert_F^2}{\det \mathbf{J}_t}\]
而在计算优化的方向时，大部分的算法都使用了局部二次近似这样的技巧。具体来说，目标函数\(E(\mathbf{x})\)可以使用泰勒展开来得到二阶近似：
\[\tilde{E}(\mathbf{x}) = E(\mathbf{x}_i) + (\mathbf{x} - \mathbf{x}_i)^T \nabla E + \frac{1}{2} (\mathbf{x} - \mathbf{x}_i)^T H_i (\mathbf{x} - \mathbf{x}_i)\]这样近似后的目标函数就是一个凸函数，可以使用各种通用求解器来进行求解。需要注意的是二次项\(H_i\)不一定严格等于目标函数的Hessian矩阵，各种算法在实现时会使用不同的技巧来估计真实的Hessian矩阵。
根据优化求解器的不同，参数化方法可以大致分为一阶方法(first-order methods)、伪牛顿法(quasi-Newton methods)以及二阶方法(second-order methods)。
在求解优化问题时需要保证三角形不会出现翻转。假设三角形上三个顶点坐标分别为\(\mathbf{u}_1\)、\(\mathbf{u}_2\)、\(\mathbf{u}_3\)，每个顶点的位移分别为\(\mathbf{v}_1 \alpha\)、\(\mathbf{v}_2 \alpha\)、\(\mathbf{v}_3 \alpha\)，其中\(\alpha\)为优化的步长。当三角形的面积发生退化时有：
\[\det \begin{bmatrix} (\mathbf{u}_2 + \mathbf{v}_2 \alpha) - (\mathbf{u}_1 + \mathbf{v}_1\alpha) \\ (\mathbf{u}_3 + \mathbf{v}_3 \alpha) - (\mathbf{u}_1 + \mathbf{v}_1\alpha) \\ \end{bmatrix} = 0\]对行列式进行展开可以得到关于\(\alpha\)的二次多项式，其最小正数根即为迭代时的最大步长。因此在整个网格上最大优化步长即为所有三角形上二次多项式的最小正根。当然在进行优化时还要保证步长满足Wolfe condition。
而常用的优化终止条件与一般的优化问题基本一致。
First-Order Methods
Block Coordinate Descent
一阶方法只利用了目标函数的一阶导数，相当于把二阶项\(H_i\)直接取0。其中分块坐标下降法(block coordinate descent, BCD)通过对优化变量进行分块的方式每一次迭代只更新一部分变量，而对于参数化的问题则可以把网格上每一个顶点作为一个block，每次只更新一个顶点坐标。
在更新顶点坐标时根据是否优化至当前最优值BCD还可以分为exact和inexact两种方法。在exact方法中每次更新顶点都会不断求解直到收敛，而inexact则只需要保证目标函数能够减少即可。实践结果显示inexact方法只需要更少的计算资源，而且往往能够得到更小的最终优化结果。
因此整个inexact BCD算法的流程可以总结如下：
Quadratic Proxy
另一种一阶方法是对总体目标函数进行分解：
\[f(\mathbf{x}) = h(\mathbf{x}) + g(\mathbf{x})\]其中\(h(\mathbf{x})\)为一个二次函数，\(h(\mathbf{x})=\frac{1}{2} \mathbf{x}^T H \mathbf{x}\)。接着对余项\(g(\mathbf{x})\)进行一阶展开，可以得到：
\[f(\mathbf{x} + \mathbf{p}) \approx h(\mathbf{x} + \mathbf{p}) + g(\mathbf{x}) + \nabla g(\mathbf{x}) \cdot \mathbf{p}\]此时计算优化方向\(\mathbf{p}\)等价于求解线性方程组：
\[H \mathbf{p} = -\nabla f(\mathbf{x})\]
在实际进行求解时还可以结合一些加速策略来加快收敛。
Scalable Locally Injective Mappings
SLIM同样是对原始目标函数进行一个凸的近似。但SLIM没有使用泰勒展开，而是使用了加权的误差来进行近似。除此之外，SLIM还要求近似后目标函数的梯度与原始目标函数的梯度相同。
整个SLIM算法的流程如下：
Quasi-Newton Methods
L-BFGS
伪牛顿法中最常见的方法是L-BFGS。首先对目标函数在\(\mathbf{x}_{n+1}\)处进行二阶泰勒展开：
\[f(\mathbf{x}) \approx f(\mathbf{x}_{n+1}) + \nabla f(\mathbf{x}_{n+1})^T (\mathbf{x} - \mathbf{x}_{n+1}) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_{n+1})^T H_{n+1} (\mathbf{x} - \mathbf{x}_{n+1})\]对\(f(\mathbf{x})\)求导后再令\(\mathbf{x} = \mathbf{x}_{n}\)可以得到：
\[\nabla f(\mathbf{x}_{n}) = \nabla f(\mathbf{x}_{n+1}) + H_{n+1} (\mathbf{x}_n - \mathbf{x}_{n+1})\]整理后得到：
\[H_{n+1} (\mathbf{x}_{n+1} - \mathbf{x}_n) = \nabla f(\mathbf{x}_{n+1}) - \nabla f(\mathbf{x}_{n})\]上式称为割线方程。记位置差异为\(\mathbf{s}_n = \mathbf{x}_{n+1} - \mathbf{x}_n\)，梯度差异为\(\mathbf{y}_n = \nabla f(\mathbf{x}_{n+1}) - \nabla f(\mathbf{x}_{n})\)，可以得到线性方程：
\[H_{n+1} \mathbf{s}_n = \mathbf{y}_n \Rightarrow \mathbf{s}_n = D_{n+1} \mathbf{y}_n\]其中\(D_{n+1} = H_{n+1}^{-1}\)。在L-BFGS中使用rank-two修正来估计\(D_{n+1}\)：
\[D_{n+1} = \bigg( I - \frac{\mathbf{y}_n \mathbf{s}_n^T}{\mathbf{s}_n^T \mathbf{y}_n} \bigg)^T D_n \bigg( I - \frac{\mathbf{y}_n \mathbf{s}_n^T}{\mathbf{s}_n^T \mathbf{y}_n} \bigg) + \frac{\mathbf{s}_n \mathbf{s}_n^T}{\mathbf{s}_n^T \mathbf{y}_n}\]更新\(D_{n+1}\)后即可利用梯度来计算优化方向：
\[\mathbf{d}_{n+1} = - D_{n+1} \nabla f(\mathbf{x}_{n+1})\]
BCQN
BCQN是对L-BFGS的改进。当几何扭曲比较大时标准L-BFGS容易产生数值不稳定的问题，因此BCQN指出可以使用网格的Laplace矩阵\(\mathbf{L}\)来计算\(\mathbf{y}_n\)：
\[\mathbf{y}_n = \mathbf{L} \mathbf{s}_n = \mathbf{L} (\mathbf{x}_{n+1} - \mathbf{x}_n)\]再结合线性加权来获得超线性的收敛：
\[\mathbf{y}_n \leftarrow (1 - \beta_n) \mathbf{y}_n + \beta_n \mathbf{L} \mathbf{s}_n\]
Second-Order Methods
二阶方法一般拥有更高的收敛速度，但在每一步的计算中往往需要更多的资源来计算Hessian矩阵以及求解线性方程组。而为了保证Hessian矩阵的正定性，还会对估计的Hessian矩阵进行一些修正。
由于大多数Hessian矩阵都可以表示为每个三角形上Hessian矩阵的和，因此只要保证每个三角形上Hessian矩阵的正定性即可。
Projected Newton
最简单的修正方法是利用特征值分解，然后把矩阵的所有小于0的特征值都设为0。不过需要注意的是矩阵的特征值分解往往比较慢，当问题的规模很大时并不推荐这种方法。
Analytic Eigen System
另一种修正方法是使用Jacobian矩阵的不变量。利用这些不变量可以解析地将Hessian矩阵投影到(半)正定矩阵空间中。
Composite Majorization
我们可以把能量函数\(f\)分解为\(h\)和\(g\)两个函数的复合，其中\(g\)将顶点坐标映射到某个空间然后在这个空间上通过\(h\)计算得到实际的能量。比如说我们可以取\(g\)为计算Jacobian矩阵特征值的过程，此时\(h\)即为利用特征值来构造所需的能量项。
显然大多数情况下\(h\)和\(g\)两个函数都是比较复杂的，不过我们可以对它们进行凹凸分解，即对一个坐标是凸函数而对另一个坐标是凹函数。
实际求解优化问题时可以利用MM(majorization minimization)将原始优化函数转换为它的一个凸上界称为majorizer，通过优化这个majorizer来实现优化。再结合凹凸分解来构造出所需的majorizer并且进行优化即可。
实际上构造出的majorizer不仅在\(x_0\)处与原函数有相同的函数值，而且还具有相同的一阶导数。
Progressive Parameterization
前面介绍过的方法其核心思想在于构造优化目标函数以及使用合适的求解器进行求解。然而大多数情况下目标函数是非常难以进行优化的，即使使用了二阶近似等技巧也需要相当多的迭代次数才会收敛。
试验结果显示，当参考三角形与初始化三角形之间存在较大的扭曲时就会导致优化比较困难。而如果二者之间的扭曲比较小，只需要少量的迭代就可以收敛。
因此渐进参数化(progressive parameterization, PP)的思想在于重新计算一个新的参考三角形，使得它与参数化之间的扭曲是有界的。然后再令构造的临时参考三角形逐步逼近原始网格，通过不断地优化和迭代就可以让最终的参数化结果逼近原始网格。由于每一步的优化都是在一个低扭曲的三角形上，整个参数化过程会比直接进行优化快很多。
