经典论文推导: As-Rigid-As-Possible(ARAP) Surface Modeling

论文As-Rigid-As-Possible Surface Modeling 发表于SGP 2007，是变形领域的经典论文，目前引用已经超过1000次。网格变形要求产生视觉合理并且大致满足物理规律的变形效果，而模型细节的保持很大程度地满足了这种需求。刚性作为一个重要的属性在保细节方面具有明显的优势，从某种意义上讲，保持模型的刚性很大程度地促进了模型表面细节的保持。

先来看看文中的变形效果：

本文主要讲一下论文中的旋转矩阵$R_i$的求解方法。

一、算法简介

如下图所示，ARAP变形算法首先定义网格顶点与1-邻域的边构成刚性变形单元，所有点的变形单元重叠地覆盖网格表面。

变形过程中假设变形单元仅仅发生旋转变换，形式化的表示如公式（1）所示，变形单元的刚性变换使得网格顶点相对于1-邻域顶点的位置保持不变，从而有效地保持了模型局部的细节。
$$p_i^{\prime }-p_j^{\prime }=R_i(p_i-p_j),\forall j\in N(i)\text{\hspace{1em}(1)}$$
其中，$R_i$ 是该变形单元对应的旋转矩阵，$p_i,p_j$ 和 $p_i^{\prime },p_j^{\prime }$ 是变形前后的网格顶点，这里点的表示是列向量。

每个变形单元的能量函数定义如下所示，通过最小化该能量函数实现模型的尽可能刚性变形：
$$E(C_i,C_i^{\prime })=\sum _{j\in N(i)}{w}_{ij}{\left|\left|(p_i^{\prime }-p_j^{\prime })-R_i(p_i-p_j)\right|\right|}^2\text{\hspace{1em}(2)}$$

ARAP变形算法总的能量函数定义为：
$$E=\sum _{i=1}^n\sum _{j\in N(i)}{w}_{ij}{\left|\left|(p_i^{\prime }-p_j^{\prime })-R_i(p_i-p_j)\right|\right|}^2\text{\hspace{1em}(3)}$$

其中，$C_i$和$C_i^{\prime }$分别表示变形前后模型顶点 $p_i$ 和 $p_i^{\prime }$ 对应的变形单元，$N(i)$ 表示 $p_i$ 的1-邻域点的索引，而 $p_j$ 和 $p_j^{\prime }$ 分别表示 $p_i$ 和 $p_i^{\prime }$ 的1-邻域顶点，$R_i$表示$C_i$到$C_i^{\prime }$的最优旋转矩阵，$w_{ij}$是边$e_{ij}=(p_i,p_j)$的权重, 定义如下所示。其中，${\alpha }_{ij}$ 和 ${\beta }_{ij}$ 如上图所示。
$$w_{i,j}=\frac{1}{2}(cot{\alpha }_{ij}+cot{\beta }_{ij})\text{\hspace{1em}(4)}$$

对于公式（3）所示的二次能量函数最小化中的每个变形单元，共有两个未知数：旋转矩阵$R$ 和变形后的点坐标 $P^{\prime }$(公式(2)中的$p_i^{\prime }$ 和 $p_j^{\prime }$). 论文采用迭代的方式进行优化求解：

• 1）固定$P^{\prime }$ 求解使得(3)最小的每个变形单元的最优旋转矩阵$R_i$。
• 2）在旋转矩阵$R_i$已知的情况下，求解使得(3)最小的变形后的顶点坐标$P^{\prime }$。
• 3）返回1），直到能量误差小于用户指定的阈值为止。

接下来，重点讲一下如何求解$R$ 和 $P^{\prime }$.

二、已知$P^{\prime }$求解旋转矩阵$R$

针对每个变形单元单独求解，令$e_{ij}^{\prime }=p_i^{\prime }-p_j^{\prime },e_{ij}=p_i-p_j$, 公式(2)可以改写为：
$$E(C_i,C_i^{\prime })=\sum _{j\in N(i)}{w}_{ij}(e_{ij}^{\prime }-R_i{e}_{ij}{)}^2=\sum _{j\in N(i)}{w}_{ij}(e_{ij}^{\prime T}{e}_{ij}^{\prime }-e_{ij}^{\prime T}{R}_i{e}_{ij}+e_{ij}^T{R}_i^T{R}_i{e}_{ij})\text{\hspace{1em}(5)}$$
因为$R$是旋转矩阵(正交矩阵)，$R_i^T{R}_i=I$, 所以上式等于：
$$E(C_i,C_i^{\prime })=\sum _{j\in N(i)}{w}_{ij}(e_{ij}^{\prime T}{e}_{ij}^{\prime }-e_{ij}^{\prime T}{R}_i{e}_{ij}+e_{ij}^T{e}_{ij})\text{\hspace{1em}(6)}$$
上式的最小值跟不含$R_i$的项无关，因此：
$$R_i=\underset{R_i}{\mathrm{argmin}}\sum _{j\in N(i)}{w}_{ij}(-e_{ij}^{\prime T}{R}_i{e}_{ij})=\underset{R_i}{\mathrm{argmax}}\sum _{j\in N(i)}{w}_{ij}{e}_{ij}^{\prime T}{R}_i{e}_{ij}\text{\hspace{1em}(7)}$$
上式可以写成矩阵的迹的形式(参见末尾示意图)进行改写：
$$R_i=\underset{R_i}{\mathrm{argmax}}Tr(\sum _{j\in N(i)}{w}_{ij}{R}_i{e}_{ij}{e}_{ij}^{\prime T})\text{\hspace{1em}(8)}$$

并将$R_i$提到外面：
$$R_i=\underset{R_i}{\mathrm{argmax}}Tr(R_i\sum _{j\in N(i)}{w}_{ij}{e}_{ij}{e}_{ij}^{\prime T})\text{\hspace{1em}(9)}$$

接下来就是如何求解这个旋转矩阵，假设$S_i={\sum }_{j\in N(i)}{w}_{ij}{e}_{ij}{e}_{ij}^{\prime T}$， 一般而言，$S_i$ 可以通过奇异值分解得到:
$$S_i=U_i{\Sigma }_i{V}_i^T\text{\hspace{1em}(10)}$$
那么(9)可以进一步写为:
$$R_i=\underset{R_i}{\mathrm{argmax}}Tr(R_i{U}_i{\Sigma }_i{V}_i^T)\text{\hspace{1em}(11)}$$
根据$Tr(AB)=Tr(BA)$，上式进一步改写：
$$R_i=\underset{R_i}{\mathrm{argmax}}Tr(V_i^T{R}_i{U}_i{\Sigma }_i)\text{\hspace{1em}(12)}$$

令：$X=V_i^T{R}_i{U}_i$, 因为此三个矩阵都是正交矩阵，所以$X$也是正交矩阵，上式可以写成：
$$R_i=\underset{R_i}{\mathrm{argmax}}Tr(X{\Sigma }_i)\text{\hspace{1em}(13)}$$
而:
$$Tr(X{\Sigma }_i)=X_{(1,1)}{\Sigma }_i(1,1)+X_{(2,2)}{\Sigma }_i(2,2)+X_{(3,3)}{\Sigma }_i(3,3)\text{\hspace{1em}(14)}$$
对于上式有两点需要注意：

1）${\Sigma }_i$保存了矩阵$S_i$的所有奇异值，所以${\Sigma }_i$的对角线上的元素均大于等于0（奇异值分解的性质决定）。
2）另外，正交矩阵中的各元素的取值$X_{(i,j)}\in [0,1]$。
因此，当X是单位阵时上式取得最大值，即${\Sigma }_i$对角线元素之和。因此只需要$X=V_i^T{R}_i{U}_i=I$ 即可，所以：
$$R_i=V_i{U}_i^T\text{\hspace{1em}(15)}$$
这就是论文给出的最终的旋转矩阵的求解方法。

三、已知$R$求解变形后的顶点$P^{\prime }$

这一步相对比较简单，直接对能量函数求导即可，具体推导就把论文中搬出来即可：
$$\begin{array}{rl}\frac{\partial E\left({\mathcal{S}}^{\prime }\right)}{\partial {\mathbf{p}}_i^{\prime }}& =\frac{\partial }{\partial {\mathbf{p}}_i^{\prime }}(\sum _{j\in \mathcal{N}(i)}{w}_{ij}{\Vert \left({\mathbf{p}}_i^{\prime }-{\mathbf{p}}_j^{\prime }\right)-{\mathbf{R}}_i\left({\mathbf{p}}_i-{\mathbf{p}}_j\right)\Vert }^2\\ & +\sum _{j\in \mathcal{N}(i)}{w}_{ji}{\Vert \left({\mathbf{p}}_j^{\prime }-{\mathbf{p}}_i^{\prime }\right)-{\mathbf{R}}_j\left({\mathbf{p}}_j-{\mathbf{p}}_i\right)\Vert }^2)=\\ & =\sum _{j\in \mathcal{N}(i)}2w_{ij}\left(\left({\mathbf{p}}_i^{\prime }-{\mathbf{p}}_j^{\prime }\right)-{\mathbf{R}}_i\left({\mathbf{p}}_i-{\mathbf{p}}_j\right)\right)+\\ & +\sum _{j\in \mathcal{N}(i)}-2w_{ji}\left(\left({\mathbf{p}}_j^{\prime }-{\mathbf{p}}_i^{\prime }\right)-{\mathbf{R}}_j\left({\mathbf{p}}_j-{\mathbf{p}}_i\right)\right).\end{array}\text{\hspace{1em}(16)}$$
因为$w_{ij}=w_{ji}$, 上式进一步改写为：
$$\frac{\partial E\left({\mathcal{S}}^{\prime }\right)}{\partial {\mathbf{p}}_i^{\prime }}=\sum _{j\in \mathcal{N}(i)}4w_{ij}\left(\left({\mathbf{p}}_i^{\prime }-{\mathbf{p}}_j^{\prime }\right)-\frac{1}{2}\left({\mathbf{R}}_i+{\mathbf{R}}_j\right)\left({\mathbf{p}}_i-{\mathbf{p}}_j\right)\right)$$
令导数为0，则可以得到：
$$\sum _{j\in \mathcal{N}(i)}{w}_{ij}\left({\mathbf{p}}_i^{\prime }-{\mathbf{p}}_j^{\prime }\right)=\sum _{j\in \mathcal{N}(i)}\frac{w_{ij}}{2}\left({\mathbf{R}}_i+{\mathbf{R}}_j\right)\left({\mathbf{p}}_i-{\mathbf{p}}_j\right)\text{\hspace{1em}(17)}$$
最终转换为一个线性方程组求解即可。

四、补充：(7)式到(8)式，矩阵到迹的表示

可令(7)式中的 $A=w_{ij}{e}_{ij}^{\prime T},B=R_i{e}_{ij}$, 其中A是一行3列，B是3行一列，对于这种行矩阵跟列矩阵的乘法可以表示为矩阵的迹，下面展示了更一般的情况: $A_{1\times n}{B}_{n\times 1}=Tr(B_{n\times 1}{A}_{1\times n})$。

参考资料:

[1] 知乎：用数学编辑3D模型（二）- As-Rigid-As-Possible
[2] 知乎：奇异值分解（SVD）
[3] 矩阵迹的性质和运算
[4] http://blog.csdn.net/seamanj/article/details/50597420 (Tr(M) >= Tr(RM))