Local/Global 参数化算法实现

关于LSCM算法实现完成的是对每个三角形尽量保证保角，即相似/共形变换，而ARAP算法的目的就是在此基础上尽量保证变换后的三角形面积相等

基本思路

假设对于空间上每个三角形 $t$t 建立局部坐标系后的坐标为 $x_t^0$xt0​ , $x_t^1$xt1​ , $x_t^2$xt2​，参数化后的坐标为 $u_t^0$ut0​ , $u_t^1$ut1​ , $u_t^2$ut2​ ,其 $jacobi$jacobi 矩阵为 $J_t$Jt​ ,

对于每个 $J_t$Jt​ 都可以进行 $SVD$SVD 分解，写成
$J=U\Sigma V^T$J=UΣVT
其中 $U$U , $V$V 为正交矩阵，$\Sigma$Σ 为对角矩阵，可以写为
$\Sigma=\left(\begin{matrix} \sigma_1&0\\ 0&\sigma_2 \end{matrix}\right)$Σ=(σ1​0​0σ2​​)
最佳的旋转变换矩阵可以记为 $L=UV^T$L=UVT 。

定义能量函数
$E(u,L)=\sum_{t=1}^TA_t\left|\left|J_t(u)-L_t\right|\right|^2$E(u,L)=t=1∑T​At​∣∣Jt​(u)−Lt​∣∣2
目标为此能量函数的最小值

而此能量函数有可写为
$E(u,L)=\frac{1}{2}\sum_{t=1}^T\sum_{i=0}^2cot(\theta_t^i)\left|\left|(u_t^i-u_t^{i+1})-L_t(x_t^i-x_t^{i+1})\right|\right|^2$E(u,L)=21​t=1∑T​i=0∑2​cot(θti​)∣∣​∣∣​(uti​−uti+1​)−Lt​(xti​−xti+1​)∣∣​∣∣​2
其中 $\theta_t^i$θti​ 为边 $(x_t^i,x_t^{i+1})$(xti​,xti+1​) 的对角，显然上述公式的 $L_t$Lt​ 也与参数化坐标 $u_t^0$ut0​ , $u_t^1$ut1​ , $u_t^2$ut2​ ,有关，这样导致能量函数的很难求解，基本思路分为两个步骤

步骤一(local phase)：固定所有的 $u$u ,求出所有对应的最佳旋转矩阵 $L_t$Lt​
步骤二(global phase)：固定所有的旋转矩阵 $L_t$Lt​ ,解上述的最小二乘系统得到对应的最佳 $u$u

重复迭代上述两个步骤即可

下面对两个步骤分别说明

local phase：

现在有着所有的三角形的局部坐标，即其参数化后的坐标，求出其 $Jacobi$Jacobi 矩阵并不是难事，利用
$\left(\begin{matrix} u_t^1-u_t^0&u_t^2-u_t^0&0\\ v_t^1-v_t^0&v_t^2-v_t^0&0\\ 1&1&1 \end{matrix}\right)= J_t\left(\begin{matrix} x_t^1-x_t^0&x_t^2-x_t^0&0\\ y_t^1-y_t^0&y_t^2-y_t^0&0\\ 1&1&1 \end{matrix}\right)$⎝⎛​ut1​−ut0​vt1​−vt0​1​ut2​−ut0​vt2​−vt0​1​001​⎠⎞​=Jt​⎝⎛​xt1​−xt0​yt1​−yt0​1​xt2​−xt0​yt2​−yt0​1​001​⎠⎞​
对于每个三角形来说其上式中 $Jacobi$Jacobi 右边的矩阵显然是满秩的，求出其逆矩阵即可求出其 $Jacobi$Jacobi 矩阵，根据上述的做法求出最佳旋转矩阵即可。

论文中给出了另一种做法

假设对于每个三角形有一个最佳的旋转矩阵使得下面的能量函数最小
$E(\pmb x,\pmb u)=\sum_{j\in\boldsymbol N(i)}w_{ij}\left|\left| (\pmb u_i-\pmb u_j)-\pmb R_i(\pmb x_i-\pmb x_j)\right|\right|^2$E(xxx,uuu)=j∈N(i)∑​wij​∣∣(uuui​−uuuj​)−RRRi​(xxxi​−xxxj​)∣∣2
定义 $\pmb e_{ij}=\pmb x_i-\pmb x_j$eeeij​=xxxi​−xxxj​ ， $\pmb e_{ij}'=\pmb u_i-\pmb u_j$eeeij′​=uuui​−uuuj​ ，上式可写为
$E(\pmb x,\pmb u)=\sum_j w_{ij}(\pmb e_{ij}’-\pmb R_i\pmb e_{ij})^T(\pmb e_{ij}’-\pmb R_i\pmb e_{ij})\\ =\sum_jw_{ij}(\pmb e_{ij}’^T\pmb e_{ij}’-2\pmb e_{ij}’^T\pmb R_i\pmb e_{ij}+\pmb e_{ij}^T\pmb e_{ij})$E(xxx,uuu)=j∑​wij​(eeeij​’−RRRi​eeeij​)T(eeeij​’−RRRi​eeeij​)=j∑​wij​(eeeij​’Teeeij​’−2eeeij​’TRRRi​eeeij​+eeeijT​eeeij​)
于是目标变为求使下面函数的最大值的 $\pmb R_i$RRRi​
${argmax}_{R_i}\sum_iw_{ij}\pmb e_{ij}’^T\pmb R_i\pmb e_{ij}\\ ={argmax}_{R_i}Tr\left(\sum_iw_{ij}\pmb R_i\pmb e_{ij}\pmb e_{ij}’^T\right)\\ ={argmax}_{R_i}Tr\left(\pmb R_i\sum_iw_{ij}\pmb e_{ij}\pmb e_{ij}’^T\right)$argmaxRi​​i∑​wij​eeeij​’TRRRi​eeeij​=argmaxRi​​Tr(i∑​wij​RRRi​eeeij​eeeij​’T)=argmaxRi​​Tr(RRRi​i∑​wij​eeeij​eeeij​’T)
定义 $\pmb S_i=\sum_{j\in \boldsymbol N(i)}w_{ij}\pmb e_{ij}\pmb e_{ij}’^T$SSSi​=∑j∈N(i)​wij​eeeij​eeeij​’T ，于是目标成为求 $Tr(\pmb R_i\pmb S_i)$Tr(RRRi​SSSi​) 最大化的 $\pmb R_i$RRRi​ （对于半正定矩阵 $\pmb M$MMM ,对任意的旋转矩阵 $\pmb R$RRR ，有 $Tr(\pmb M)\geq Tr(\pmb R\pmb M)$Tr(MMM)≥Tr(RRRMMM) ）显然对 $\pmb S_i$SSSi​ 进行SVD分解，$\pmb S_i=\pmb U_i\pmb\Sigma_i\pmb V_i^T$SSSi​=UUUi​ΣΣΣi​VVViT​ ，不难知道当 $\pmb R_i=\pmb V_i\pmb U_i^T$RRRi​=VVVi​UUUiT​ 时满足条件。上述的 $w_{ij}$wij​ 为 $Laplace$Laplace 算子，一般定义为 $w_{ij}=\frac{1}{2}(cot(\theta_{ij})+cot(\theta_{ji}))$wij​=21​(cot(θij​)+cot(θji​))

事实上，通过解逆矩阵的方式也能得到不错的效果，不过效率相对低下

global phase:

上述提到能量函数公式可以进一步以三角网格的半边表示
$E(u,L)=\frac{1}{2}\sum_{t=1}^T\sum_{i=0}^2cot(\theta_t^i)\left|\left|(u_t^i-u_t^{i+1})-L_t(x_t^i-x_t^{i+1})\right|\right|^2\\ =\frac{1}{2}\sum_{(i,j)\in he}cot(\theta_{ij})\left|\left|(u_i-u_{j})-L_{t(i,j)}(x_i-x_j)\right|\right|^2$E(u,L)=21​t=1∑T​i=0∑2​cot(θti​)∣∣​∣∣​(uti​−uti+1​)−Lt​(xti​−xti+1​)∣∣​∣∣​2=21​(i,j)∈he∑​cot(θij​)∣∣​∣∣​(ui​−uj​)−Lt(i,j)​(xi​−xj​)∣∣​∣∣​2
对此能量求解最小值，如上篇所说，设定其梯度(每个导数)为0求极值即可，这样可以得到下方的线性系统
$\sum_{j\in \boldsymbol N(i)}[cot(\theta_{ij})+cot(\theta_{ji})](u_i-u_j)\\ =\sum_{j\in \boldsymbol N(i)}[cot(\theta_{ij})L_{t(i,j)}+cot(\theta_{ji})L_{t(j,i)}](x_i-x_j)\\ =\sum_{j\in \boldsymbol N(i)}cot(\theta_{ij})L_{t(i,j)}(x_i-x_j)+cot(\theta_{ji})L_{t(j,i)}(x_i-x_j) \\(i=1,...,n)$j∈N(i)∑​[cot(θij​)+cot(θji​)](ui​−uj​)=j∈N(i)∑​[cot(θij​)Lt(i,j)​+cot(θji​)Lt(j,i)​](xi​−xj​)=j∈N(i)∑​cot(θij​)Lt(i,j)​(xi​−xj​)+cot(θji​)Lt(j,i)​(xi​−xj​)(i=1,...,n)
显然此线性系统的左矩阵为固定矩阵，可提前求出存放，其主对角元素为所有半边所对应的余切值之和，第 $i$i 行 $j$j 列数值为其半边所在的边对应的两个角的余切值之和的相反数

其右边矩阵构造稍微复杂，但不难看出一条半边 $(i,j)$(i,j) 对其第 $i$i 行与第 $j$j 行的贡献为相反数，构造思路就是在遍历每一个三角形的每一条边时，计算 $cot(\theta_{ij})L_{t(i,j)}(x_i-x_j)$cot(θij​)Lt(i,j)​(xi​−xj​) 加在右矩阵的第 $i$i 行，并在第 $j$j 行加上其负值即可。

算法步骤

1. 得到一个初始的参数化坐标，最小二乘与凸组合均可
2. 保存每个三角形的边向量，及对应的余切值，提前构造global phase线性系统的左矩阵
3. local phase：根据当前的参数化坐标计算所有三角形的最佳旋转矩阵
4. global phase：根据计算得到的最佳旋转矩阵解线性系统得到新的参数化坐标
5. 迭代3，4步，直到能量相差变化不大即可

算法结果

测试模型

	初始	迭代1次	迭代10次	迭代20次	最终结果
凸组合
LSCM

可以看到的是初始的参数化坐标对local/global算法的最终结果并没有什么影响，不同的是，在本例中以凸组合为初始化的参数化坐标在经历第一次迭代后的效果相当糟糕，最终经历了100+次迭代到了最终结果，以LSCM为初始的参数化坐标在经历第一次迭代后的效果还算不错，最终经历了30+次迭代到课最终结果。

参考文献

[1] Liu L, Zhang L, Xu Y, et al. A local/global approach to mesh parameterization[C]. symposium on geometry processing, 2008, 27(5): 1495-1504.
[2] Sorkine O, Alexa M. As-rigid-as-possible surface modeling[C]. symposium on geometry processing, 2007: 109-116.
[3] Pinkall U, Polthier K. Computing discrete minimal surfaces and their conjugates[J]. Experimental Mathematics, 1993, 2(1): 15-36.