论文As-Rigid-As-Possible Surface Modeling 发表于SGP 2007，是变形领域的经典论文，目前引用已经超过1000次。网格变形要求产生视觉合理并且大致满足物理规律的变形效果，而模型细节的保持很大程度地满足了这种需求。刚性作为一个重要的属性在保细节方面具有明显的优势，从某种意义上讲，保持模型的刚性很大程度地促进了模型表面细节的保持。
先来看看文中的变形效果：
本文主要讲一下论文中的旋转矩阵 R i R_i Ri的求解方法。
一、算法简介
如下图所示，ARAP变形算法首先定义网格顶点与1-邻域的边构成刚性变形单元，所有点的变形单元重叠地覆盖网格表面。
变形过程中假设变形单元仅仅发生旋转变换，形式化的表示如公式（1）所示，变形单元的刚性变换使得网格顶点相对于1-邻域顶点的位置保持不变，从而有效地保持了模型局部的细节。
p i ′ − p j ′ = R i ( p i − p j ) , ∀ j ∈ N ( i ) (1) p'_i-p'_j= R_i(p_i-p_j), \forall j \in N(i)\tag1 pi′−pj′=Ri(pi−pj),∀j∈N(i)(1)
其中， R i R_i Ri 是该变形单元对应的旋转矩阵， p i , p j p_i,p_j pi,pj 和 p i ′ , p j ′ p'_i,p'_j pi′,pj′ 是变形前后的网格顶点，这里点的表示是列向量。
每个变形单元的能量函数定义如下所示，通过最小化该能量函数实现模型的尽可能刚性变形：
E ( C i , C i ′ ) = ∑ j ∈ N ( i ) w i j ∣ ∣ ( p i ′ − p j ′ ) − R i ( p i − p j ) ∣ ∣ 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 = 1 n ∑ j ∈ N ( i ) w i j ∣ ∣ ( p i ′ − p j ′ ) − R i ( p i − p j ) ∣ ∣ 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)
其中， C i C_i Ci和 C i ′ C_i' Ci′分别表示变形前后模型顶点 p i p_i pi 和 p i ′ p'_i pi′ 对应的变形单元， N ( i ) N(i) N(i) 表示 p i p_i pi 的1-邻域点的索引，而 p j p_j pj 和 p j ′ p'_j pj′ 分别表示 p i p_i pi 和 p i ′ p'_i pi′ 的1-邻域顶点， R i R_i Ri表示 C i C_i Ci到 C i ′ C_i' Ci′的最优旋转矩阵， w i j w_{ij} wij是边 e i j = ( p i , p j ) e_{ij}=(p_i, p_j) eij=(pi,pj)的权重, 定义如下所示。其中， α i j \alpha_{ij} αij 和 β i j \beta_{ij} βij 如上图所示。
w i , j = 1 2 ( c o t α i j + c o t β i j ) (4) w_{i,j} = \frac{1}{2}(cot\alpha_{ij} + cot \beta_{ij})\tag4 wi,j=21(cotαij+cotβij)(4)
对于公式（3）所示的二次能量函数最小化中的每个变形单元，共有两个未知数：旋转矩阵 R R R 和变形后的点坐标 P ′ P' P′(公式(2)中的 p i ′ p'_i pi′ 和 p j ′ p'_j pj′). 论文采用迭代的方式进行优化求解：
- 1）固定 P ′ P' P′ 求解使得(3)最小的每个变形单元的最优旋转矩阵 R i R_i Ri。
- 2）在旋转矩阵 R i R_i Ri已知的情况下，求解使得(3)最小的变形后的顶点坐标 P ′ P' P′。
- 3）返回1），直到能量误差小于用户指定的阈值为止。
接下来，重点讲一下如何求解 R R R 和 P ′ P' P′.
二、已知 P ′ P' P′求解旋转矩阵 R R R
针对每个变形单元单独求解，令 e i j ′ = p i ′ − p j ′ , e i j = p i − p j e'_{ij} = p'_i-p'_j, e_{ij} = p_i-p_j eij′=pi′−pj′,eij=pi−pj, 公式(2)可以改写为：
E ( C i , C i ′ ) = ∑ j ∈ N ( i ) w i j ( e i j ′ − R i e i j ) 2 = ∑ j ∈ N ( i ) w i j ( e i j ′ T e i j ′ − e i j ′ T R i e i j + e i j T R i T R i e i j ) (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})\tag5 E(Ci,Ci′)=j∈N(i)∑wij(eij′−Rieij)2=j∈N(i)∑wij(eij′Teij′−eij′TRieij+eijTRiTRieij)(5)
因为 R R R是旋转矩阵(正交矩阵)， R i T R i = I R^T_iR_i = I RiTRi=I, 所以上式等于：
E ( C i , C i ′ ) = ∑ j ∈ N ( i ) w i j ( e i j ′ T e i j ′ − e i j ′ T R i e i j + e i j T e i j ) (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})\tag6 E(Ci,Ci′)=j∈N(i)∑wij(eij′Teij′−eij′TRieij+eijTeij)(6)
上式的最小值跟不含 R i R_i Ri的项无关，因此：
R i = arg min R i ∑ j ∈ N ( i ) w i j ( − e i j ′ T R i e i j ) = arg max R i ∑ j ∈ N ( i ) w i j e i j ′ T R i e i j (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}\tag7 Ri=Riargminj∈N(i)∑wij(−eij′TRieij)=Riargmaxj∈N(i)∑wijeij′TRieij(7)
上式可以写成矩阵的迹的形式(参见末尾示意图)进行改写：
R i = arg max R i T r ( ∑ j ∈ N ( i ) w i j R i e i j e i j ′ T ) (8) R_i =\argmax \limits_{R_i}Tr(\sum_{j\in N(i)}w_{ij}R_ie_{ij}e'^T_{ij}) \tag8 Ri=RiargmaxTr(j∈N(i)∑wijRieijeij′T)(8)
并将 R i R_i Ri提到外面：
R i = arg max R i T r ( R i ∑ j ∈ N ( i ) w i j e i j e i j ′ T ) (9) R_i =\argmax \limits_{R_i}Tr(R_i\sum_{j\in N(i)}w_{ij}e_{ij}e'^T_{ij})\tag9 Ri=RiargmaxTr(Rij∈N(i)∑wijeijeij′T)(9)
接下来就是如何求解这个旋转矩阵，假设 S i = ∑ j ∈ N ( i ) w i j e i j e i j ′ T S_i = \sum_{j\in N(i)}w_{ij}e_{ij}e'^T_{ij} Si=∑j∈N(i)wijeijeij′T， 一般而言， S i S_i Si 可以通过奇异值分解得到:
S i = U i Σ i V i T (10) S_i = U_i\Sigma_iV^T_i\tag{10} Si=UiΣiViT(10)
那么(9)可以进一步写为:
R i = arg max R i T r ( R i U i Σ i V i T ) (11) R_i =\argmax \limits_{R_i}Tr(R_iU_i\Sigma_iV^T_i) \tag{11} Ri=RiargmaxTr(RiUiΣiViT)(11)
根据 T r ( A B ) = T r ( B A ) Tr(AB) = Tr(BA) Tr(AB)=Tr(BA)，上式进一步改写：
R i = arg max R i T r ( V i T R i U i Σ i ) (12) R_i =\argmax \limits_{R_i}Tr(V^T_iR_iU_i\Sigma_i) \tag{12} Ri=RiargmaxTr(ViTRiUiΣi)(12)
令： X = V i T R i U i X = V^T_iR_iU_i X=ViTRiUi, 因为此三个矩阵都是正交矩阵，所以 X X X也是正交矩阵，上式可以写成：
R i = arg max R i T r ( X Σ i ) (13) R_i =\argmax \limits_{R_i}Tr(X\Sigma_i) \tag{13} Ri=RiargmaxTr(XΣi)(13)
而:
T r ( 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保存了矩阵 S i S_i Si的所有奇异值，所以 Σ 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 = V i T R i U i = I X= V^T_iR_iU_i=I X=ViTRiUi=I 即可，所以：
R i = V i U i T (15) R_i = V_iU^T_i\tag{15} Ri=ViUiT(15)
这就是论文给出的最终的旋转矩阵的求解方法。
三、已知 R R R求解变形后的顶点 P ′ P' P′
这一步相对比较简单，直接对能量函数求导即可，具体推导就把论文中搬出来即可：
∂ E ( S ′ ) ∂ p i ′ = ∂ ∂ p i ′ ( ∑ j ∈ N ( i ) w i j ∥ ( p i ′ − p j ′ ) − R i ( p i − p j ) ∥ 2 + ∑ j ∈ N ( i ) w j i ∥ ( p j ′ − p i ′ ) − R j ( p j − p i ) ∥ 2 ) = = ∑ j ∈ N ( i ) 2 w i j ( ( p i ′ − p j ′ ) − R i ( p i − p j ) ) + + ∑ j ∈ N ( i ) − 2 w j i ( ( p j ′ − p i ′ ) − R j ( p j − p i ) ) . (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)
因为 w i j = w j i w_{ij} = w_{ji} wij=wji, 上式进一步改写为：
∂ E ( S ′ ) ∂ p i ′ = ∑ j ∈ N ( i ) 4 w i j ( ( p i ′ − p j ′ ) − 1 2 ( R i + R j ) ( p i − p j ) ) \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 ) w i j ( p i ′ − p j ′ ) = ∑ j ∈ N ( i ) w i j 2 ( R i + R j ) ( p i − p j ) (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 = w i j e i j ′ T , B = R i e i j A = w_{ij}e'^T_{ij}, B = R_ie_{ij} A=wijeij′T,B=Rieij, 其中A是一行3列，B是3行一列，对于这种行矩阵跟列矩阵的乘法可以表示为矩阵的迹，下面展示了更一般的情况: A 1 × n B n × 1 = T r ( B n × 1 A 1 × 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))
