目录

刚性变换的问题描述

最优平移向量求解

最优旋转矩阵求解

反射矩阵消除

基于SVD刚性变换矩阵计算流程总结

---

刚性变换的问题描述

令P={p_1,p_2,...,p_n}和Q={q_1,q_2,...,q_n}是R^d空间内的两组对应的点。希望找到一个刚性的变换，在最小二乘的意义上最优地对齐两个点集，也就是说，寻求一个旋转矩阵R和一个平移向量t来满足： (R,t)=argmin┬RϵSO(d),tϵR^d∑_i=1^n▒ω_i‖(Rp_i+t)−q_i‖^2                （1） 其中ω_i>0是每个对点的权重。

最优平移向量求解

假设R被固定并定义F(t)=∑_i=1^n▒ω_i‖(Rp_i+t)−q_i‖^2。参考The Matrix Cookbook可知∂‖x‖_2/∂x=∂(x^Tx)/∂x=2x。可以通过取F(t)相对于t的导数并计算它的根来找到最优平移向量t： 0=∂F/∂x=∑_i=1^n▒2ω_i(Rp_i+t−q_i)                                        （2） =2t(∑_i=1^n▒ω_i)+2R(∑_i=1^n▒ω_ip_i)−2∑_i=1^n▒ω_iq_i （3） 定义 p ̅=∑_i=1^n▒ω_ip_i/∑_i=1^n▒ω_i，q ̅=∑_i=1^n▒ω_iq_i/∑_i=1^n▒ω_i                                            （4） 通过重新排列上面公式的项，可以得到： t=q ̅−Rp ̅                                                          （5） 换句话说，最优平移向量t将被变换的点集P的加权质心映射到点集Q的加权质心。将最优平移向量t代入目标函数： ∑_i=1^n▒ω_i‖(Rp_i+t)−q_i‖^2=∑_i=1^n▒ω_i‖Rp_i+q ̅−Rp ̅−q_i‖^2 （6）                                     =∑_i=1^n▒ω_i‖R(p_i−q_i)−((q_i) ̅−q ̅)‖^2                        （7）  所以寻找的最佳旋转矩阵R，使其满足： R= argmin┬RϵSO(d),tϵR^d∑_i=1^n▒ω_i‖Rx_i−y_i‖^2 （8）

最优旋转矩阵求解

对公式（8）重新整理可得： ‖Rx_i−y_i‖^2=(Rx_i−y_i)^T(Rx_i−y_i)=(x_i^TR^T−y_i^T)(Rx_i−y_i) =x_i^TR^TRx_i−y_i^T Rx_i−x_i^TR^Ty_i+y_i^Ty_i =x_i^Tx_i−y_i^T Rx_i−x_i^TR^Ty_i+y_i^Ty_i (9) 通过旋转矩阵的正交性RR^T=R^TR=Ⅰ（Ⅰ是单位矩阵）。 需要注意的是，x_i^TR^Ty_i是一个标量：x_i^T维度为1×d，R^T维度为d×d并且y_i的维度为d×1。对于任何标量a，通常都有a=a^T，因此 x_i^TR^Ty_i=(x_i^TR^Ty_i)^T=y_i^TRx_i (10) 因此有 ‖Rx_i−y_i‖^2=x_i^Tx_i−2y_i^T Rx_i+y_i^Ty_i (11) 将公式（11）替换公式（8）： argmin┬RϵSO(d)∑_i=1^n▒ω_i‖Rx_i−y_i‖^2=argmin┬RϵSO(d)∑_i=1^n▒ω_i(x_i^Tx_i−2y_i^T Rx_i+y_i^Ty_i) =argmin┬RϵSO(d)(∑_i^n▒ω_ix_i^Tx_i−2∑_i^n▒ω_iy_i^T Rx_i+∑_i^n▒ω_iy_i^Ty_i) =argmin┬RϵSO(d)(−2∑_i^n▒ω_iy_i^T Rx_i) (12)

最后一步（移除∑_i=1^n▒ω_ix_i^Tx_i和∑_i=1^n▒ω_iy_i^Ty_i）成立，因为这些表达式不依赖于R，所以移除他们不会影响最小值。最小化表达式乘以标量也是如此，所以有 argmin┬RϵSO(d)(−2∑_i^n▒ω_iy_i^T Rx_i)=argmax┬RϵSO(d)∑_i=1^n▒ω_iy_i^T Rx_i (13) 将上式可以改为矩阵形式如下： ∑_i=1^n▒ω_iy_i^T Rx_i=tr(WY^TRX) (14) 其中W=diag(ω_1,ω_2,…,ω_n)是带加权对角元素ω_i的n×n对角矩阵；Y是一个以y_i为列的d×n矩阵，X是一个以x_i为列的d×n矩阵。在这里提醒读者，方阵的迹是对角线上元素的和：tr(A)=∑_i=1^n▒a_ii。下图给出了代数操作的说明。

反射矩阵消除

上面得到的矩阵并不能保证是一个旋转矩阵，可能为反射矩阵（Reflection matrix），此时可以通过验证VU^T的行列式来判断到底是旋转（det(VU^T)=1）还是反射（det(VU^T)=−1）。为了确保是旋转矩阵，这时需要对公式（21）进一步处理。 假设det(VU^T)=−1，则限制R为旋转矩阵就意味着M=V^TRU为反射矩阵，于是找到一个反射矩阵M来最大化： tr(∑M)=σ_1m_11+σ_2m_22+…+σ_dm_dd=: f(m_11,m_22,…,m_dd) (22) 其中f是以m_11,m_22,…,m_dd为变量的线性函数，由于m_ii∈[−1,1]，其极大值在其定义域的边界处。于是当∀i,m_ii=1时，f取得最大值，但是此时的R为反射矩阵，所以并不能这样取值。然后来看第二个极大值点(1,1,…,−1)，有： f=tr(∑M)=σ_1+σ_2+…+σ_d−1−σ_d (23) 这个值大于任何其他的自变量取值(±1,±1,…,±1)除了(1,1,…,1)的组合，因为奇异值是经过排序的，σ_d是最小的奇异值。

基于SVD刚性变换矩阵计算流程总结