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

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

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

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

一、算法简介

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

变形过程中假设变形单元仅仅发生旋转变换，形式化的表示如公式（1）所示，变形单元的刚性变换使得网格顶点相对于1-邻域顶点的位置保持不变，从而有效地保持了模型局部的细节。
pi′−pj′=Ri(pi−pj),∀j∈N(i)(1)p'_i-p'_j= R_i(p_i-p_j), \forall j \in N(i)\tag1pi′−pj′=Ri(pi−pj),∀j∈N(i)(1)
其中，RiR_iRi 是该变形单元对应的旋转矩阵，pi,pjp_i,p_jpi,pj 和 pi′,pj′p'_i,p'_jpi′,pj′ 是变形前后的网格顶点，这里点的表示是列向量。

每个变形单元的能量函数定义如下所示，通过最小化该能量函数实现模型的尽可能刚性变形：
E(Ci,Ci′)=∑j∈N(i)wij∣∣(pi′−pj′)−Ri(pi−pj)∣∣2(2)E(C_i,C_i' )= \sum_{j\in N(i)} w_{ij} \left| \left|(p'_i-p'_j) - R_i(p_i-p_j)\right|\right|^2 \tag2 E(Ci,Ci′)=j∈N(i)∑wij∣∣∣∣(pi′−pj′)−Ri(pi−pj)∣∣∣∣2(2)

ARAP变形算法总的能量函数定义为：
E=∑i=1n∑j∈N(i)wij∣∣(pi′−pj′)−Ri(pi−pj)∣∣2(3)E= \sum^n_{i=1}\sum_{j\in N(i)} w_{ij} \left| \left|(p'_i-p'_j) - R_i(p_i-p_j)\right|\right|^2 \tag3 E=i=1∑nj∈N(i)∑wij∣∣∣∣(pi′−pj′)−Ri(pi−pj)∣∣∣∣2(3)

其中，CiC_iCi和Ci′C_i'Ci′分别表示变形前后模型顶点 pip_ipi 和 pi′p'_ipi′ 对应的变形单元，N(i)N(i)N(i) 表示 pip_ipi 的1-邻域点的索引，而 pjp_jpj 和 pj′p'_jpj′ 分别表示 pip_ipi 和 pi′p'_ipi′ 的1-邻域顶点，RiR_iRi表示CiC_iCi到Ci′C_i'Ci′的最优旋转矩阵，wijw_{ij}wij是边eij=(pi,pj)e_{ij}=(p_i, p_j)eij=(pi,pj)的权重, 定义如下所示。其中，αij\alpha_{ij}αij 和 βij\beta_{ij}βij 如上图所示。
wi,j=12(cotαij+cotβij)(4)w_{i,j} = \frac{1}{2}(cot\alpha_{ij} + cot \beta_{ij})\tag4wi,j=21(cotαij+cotβij)(4)

对于公式（3）所示的二次能量函数最小化中的每个变形单元，共有两个未知数：旋转矩阵RRR 和变形后的点坐标 P′P'P′(公式(2)中的pi′p'_ipi′ 和 pj′p'_jpj′). 论文采用迭代的方式进行优化求解：

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

接下来，重点讲一下如何求解RRR 和 P′P'P′.

二、已知P′P'P′求解旋转矩阵RRR

针对每个变形单元单独求解，令eij′=pi′−pj′,eij=pi−pje'_{ij} = p'_i-p'_j, e_{ij} = p_i-p_jeij′=pi′−pj′,eij=pi−pj, 公式(2)可以改写为：
E(Ci,Ci′)=∑j∈N(i)wij(eij′−Rieij)2=∑j∈N(i)wij(eij′Teij′−eij′TRieij+eijTRiTRieij)(5)E(C_i,C_i' )=\sum_{j\in N(i)}w_{ij}(e'_{ij}-R_ie_{ij})^2 =\sum_{j\in N(i)}w_{ij}(e'^T_{ij}e'_{ij}-e'^T_{ij}R_ie_{ij}+e^T_{ij}R^T_iR_ie_{ij})\tag5E(Ci,Ci′)=j∈N(i)∑wij(eij′−Rieij)2=j∈N(i)∑wij(eij′Teij′−eij′TRieij+eijTRiTRieij)(5)
因为RRR是旋转矩阵(正交矩阵)，RiTRi=IR^T_iR_i = IRiTRi=I, 所以上式等于：
E(Ci,Ci′)=∑j∈N(i)wij(eij′Teij′−eij′TRieij+eijTeij)(6)E(C_i,C_i' )=\sum_{j\in N(i)}w_{ij}(e'^T_{ij}e'_{ij}-e'^T_{ij}R_ie_{ij}+e^T_{ij}e_{ij})\tag6E(Ci,Ci′)=j∈N(i)∑wij(eij′Teij′−eij′TRieij+eijTeij)(6)
上式的最小值跟不含RiR_iRi的项无关，因此：
Ri=arg min⁡Ri∑j∈N(i)wij(−eij′TRieij)=arg max⁡Ri∑j∈N(i)wijeij′TRieij(7)R_i = \argmin \limits_{R_i}\sum_{j\in N(i)}w_{ij}(-e'^T_{ij}R_ie_{ij}) = \argmax \limits_{R_i}\sum_{j\in N(i)}w_{ij}e'^T_{ij}R_ie_{ij}\tag7Ri=Riargminj∈N(i)∑wij(−eij′TRieij)=Riargmaxj∈N(i)∑wijeij′TRieij(7)
上式可以写成矩阵的迹的形式(参见末尾示意图)进行改写：
Ri=arg max⁡RiTr(∑j∈N(i)wijRieijeij′T)(8)R_i =\argmax \limits_{R_i}Tr(\sum_{j\in N(i)}w_{ij}R_ie_{ij}e'^T_{ij}) \tag8Ri=RiargmaxTr(j∈N(i)∑wijRieijeij′T)(8)

并将RiR_iRi提到外面：
Ri=arg max⁡RiTr(Ri∑j∈N(i)wijeijeij′T)(9)R_i =\argmax \limits_{R_i}Tr(R_i\sum_{j\in N(i)}w_{ij}e_{ij}e'^T_{ij})\tag9Ri=RiargmaxTr(Rij∈N(i)∑wijeijeij′T)(9)

接下来就是如何求解这个旋转矩阵，假设Si=∑j∈N(i)wijeijeij′TS_i = \sum_{j\in N(i)}w_{ij}e_{ij}e'^T_{ij}Si=∑j∈N(i)wijeijeij′T，一般而言，SiS_iSi 可以通过奇异值分解得到:
Si=UiΣiViT(10)S_i = U_i\Sigma_iV^T_i\tag{10}Si=UiΣiViT(10)
那么(9)可以进一步写为:
Ri=arg max⁡RiTr(RiUiΣiViT)(11)R_i =\argmax \limits_{R_i}Tr(R_iU_i\Sigma_iV^T_i) \tag{11}Ri=RiargmaxTr(RiUiΣiViT)(11)
根据Tr(AB)=Tr(BA)Tr(AB) = Tr(BA)Tr(AB)=Tr(BA)，上式进一步改写：
Ri=arg max⁡RiTr(ViTRiUiΣi)(12)R_i =\argmax \limits_{R_i}Tr(V^T_iR_iU_i\Sigma_i) \tag{12}Ri=RiargmaxTr(ViTRiUiΣi)(12)

令：X=ViTRiUiX = V^T_iR_iU_iX=ViTRiUi, 因为此三个矩阵都是正交矩阵，所以XXX也是正交矩阵，上式可以写成：
Ri=arg max⁡RiTr(XΣi)(13)R_i =\argmax \limits_{R_i}Tr(X\Sigma_i) \tag{13}Ri=RiargmaxTr(XΣi)(13)
而:
Tr(XΣi)=X(1,1)Σi(1,1)+X(2,2)Σi(2,2)+X(3,3)Σi(3,3)(14)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)\tag {14}Tr(XΣi)=X(1,1)Σi(1,1)+X(2,2)Σi(2,2)+X(3,3)Σi(3,3)(14)
对于上式有两点需要注意：

1）Σi\Sigma_iΣi保存了矩阵SiS_iSi的所有奇异值，所以Σi\Sigma_iΣi的对角线上的元素均大于等于0（奇异值分解的性质决定）。
2）另外，正交矩阵中的各元素的取值X(i,j)∈[0,1]X_{(i,j)} \in [0,1]X(i,j)∈[0,1]。
因此，当X是单位阵时上式取得最大值，即Σi\Sigma_iΣi对角线元素之和。因此只需要X=ViTRiUi=IX= V^T_iR_iU_i=IX=ViTRiUi=I 即可，所以：
Ri=ViUiT(15)R_i = V_iU^T_i\tag{15}Ri=ViUiT(15)
这就是论文给出的最终的旋转矩阵的求解方法。

三、已知RRR求解变形后的顶点P′P'P′

这一步相对比较简单，直接对能量函数求导即可，具体推导就把论文中搬出来即可：
∂E(S′)∂pi′=∂∂pi′(∑j∈N(i)wij∥(pi′−pj′)−Ri(pi−pj)∥2+∑j∈N(i)wji∥(pj′−pi′)−Rj(pj−pi)∥2)==∑j∈N(i)2wij((pi′−pj′)−Ri(pi−pj))++∑j∈N(i)−2wji((pj′−pi′)−Rj(pj−pi)).(16)\begin{aligned} \frac{\partial E\left(\mathcal{S}^{\prime}\right)}{\partial \mathbf{p}_{i}^{\prime}} &=\frac{\partial}{\partial \mathbf{p}_{i}^{\prime}}\left(\sum_{j \in \mathcal{N}(i)} w_{i j}\left\|\left(\mathbf{p}_{i}^{\prime}-\mathbf{p}_{j}^{\prime}\right)-\mathbf{R}_{i}\left(\mathbf{p}_{i}-\mathbf{p}_{j}\right)\right\|^{2}\right.\\ &\left.+\sum_{j \in \mathcal{N}(i)} w_{j i}\left\|\left(\mathbf{p}_{j}^{\prime}-\mathbf{p}_{i}^{\prime}\right)-\mathbf{R}_{j}\left(\mathbf{p}_{j}-\mathbf{p}_{i}\right)\right\|^{2}\right)=\\ &=\sum_{j \in \mathcal{N}(i)} 2 w_{i j}\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)}-2 w_{j i}\left(\left(\mathbf{p}_{j}^{\prime}-\mathbf{p}_{i}^{\prime}\right)-\mathbf{R}_{j}\left(\mathbf{p}_{j}-\mathbf{p}_{i}\right)\right) . \end{aligned}\tag{16} ∂pi′∂E(S′)=∂pi′∂⎝⎛j∈N(i)∑wij∥∥(pi′−pj′)−Ri(pi−pj)∥∥2+j∈N(i)∑wji∥∥(pj′−pi′)−Rj(pj−pi)∥∥2⎠⎞==j∈N(i)∑2wij((pi′−pj′)−Ri(pi−pj))++j∈N(i)∑−2wji((pj′−pi′)−Rj(pj−pi)).(16)
因为wij=wjiw_{ij} = w_{ji}wij=wji, 上式进一步改写为：
∂E(S′)∂pi′=∑j∈N(i)4wij((pi′−pj′)−12(Ri+Rj)(pi−pj))\frac{\partial E\left(\mathcal{S}^{\prime}\right)}{\partial \mathbf{p}_{i}^{\prime}}=\sum_{j \in \mathcal{N}(i)} 4 w_{i j}\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)∂pi′∂E(S′)=j∈N(i)∑4wij((pi′−pj′)−21(Ri+Rj)(pi−pj))
令导数为0，则可以得到：
∑j∈N(i)wij(pi′−pj′)=∑j∈N(i)wij2(Ri+Rj)(pi−pj)(17)\sum_{j \in \mathcal{N}(i)} w_{i j}\left(\mathbf{p}_{i}^{\prime}-\mathbf{p}_{j}^{\prime}\right)=\sum_{j \in \mathcal{N}(i)} \frac{w_{i j}}{2}\left(\mathbf{R}_{i}+\mathbf{R}_{j}\right)\left(\mathbf{p}_{i}-\mathbf{p}_{j}\right) \tag{17}j∈N(i)∑wij(pi′−pj′)=j∈N(i)∑2wij(Ri+Rj)(pi−pj)(17)
最终转换为一个线性方程组求解即可。

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

可令(7)式中的 A=wijeij′T,B=RieijA = w_{ij}e'^T_{ij}, B = R_ie_{ij}A=wijeij′T,B=Rieij, 其中A是一行3列，B是3行一列，对于这种行矩阵跟列矩阵的乘法可以表示为矩阵的迹，下面展示了更一般的情况: A1×nBn×1=Tr(Bn×1A1×n)A_{1\times n}B_{n\times 1} = Tr(B_{n\times 1}A_{1\times n})A1×nBn×1=Tr(Bn×1A1×n)。

参考资料:

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