YSYSIMONCODE. ART. SCIENCE. VISION.

仅包含旋转的 SVD Rotation-Variant SVD

用 SVD 和反射矩阵构造仅包含旋转的极分解。

Aug 16, 2024约 4 分钟Simulationysysimon
SimulationDeformationSVDPolar Decomposition

仅包含旋转的 SVD Rotation-Variant SVD

对于传统的 极分解 Polar decomposition F=R^S^\mathbf{F} = \mathbf{\hat{R}}\mathbf{\hat{S}}，仅要求 R^\mathbf{\hat{R}} 是酉变换，它不但包含了旋转，还有可能包含我们不想要的 反射 reflection ，所以我们必须想办法构造一种特殊的极分解，其仅包含旋转

一种直接计算矩阵 F\mathbf{F} 极分解的方法是首先计算它的 SVD

F=U^Σ^V^T\mathbf{F} = \mathbf{\hat{U}}\mathbf{\hat{\Sigma}}\mathbf{\hat{V}}^T

然后将 SVD 的结果重新组合成 R^\mathbf{\hat{R}} 和 S^\mathbf{\hat{S}}

R^=U^V^TS^=V^Σ^V^T\begin{align*} \mathbf{\hat{R}} &= \mathbf{\hat{U}}\mathbf{\hat{V}}^T \\ \mathbf{\hat{S}} &= \mathbf{\hat{V}}\mathbf{\hat{\Sigma}}\mathbf{\hat{V}}^T \end{align*}

我们的目标是找到 U^\mathbf{\hat{U}} 和 V^\mathbf{\hat{V}} 的纯旋转版本 U\mathbf{U} 和 V\mathbf{V}，然后我们可以组合它们得到最终纯旋转 R=UVT\mathbf{R} = \mathbf{U}\mathbf{V}^T

一种检测 R^\mathbf{\hat{R}} 是否含有反射的简单办法是求 R^\mathbf{\hat{R}} 的 行列式 determinant，也就是 det⁡R^\det\mathbf{\hat{R}}

  1. 不含有反射: 若 R^\mathbf{\hat{R}} 是纯旋转，则 det⁡R^=1\det\mathbf{\hat{R}} = 1 且奇异值矩阵是 Σ=[100010001]\mathbf{\Sigma} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}
  2. 含有反射: 若 R^\mathbf{\hat{R}} 含有反射，则 det⁡R^=−1\det\mathbf{\hat{R}} = -1 且奇异值矩阵必然存在一个奇异值为负，即 Σ=[10001000−1]\mathbf{\Sigma} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1 \end{bmatrix}

不存在所谓 “双重反射”，因为两次反射相当于一次旋转，如果你反射某个对象两次，这与绕着一个轴旋转 π\pi 度 (180度) 是等价的，这种旋转是一种有效的旋转操作，并不是反射，如果 F\mathbf{F} 包含这种旋转，则会在分解之后的 U^\mathbf{\hat{U}} 和 V^\mathbf{\hat{V}} 中体现，而非 Σ^\mathbf{\hat{\Sigma}} 中

有了这些信息，我们可以构造一个反射矩阵 L\mathbf{L}

L=[10001000det⁡(U^V^T)]\mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \det\left(\mathbf{\hat{U}}\mathbf{\hat{V}}^T \right) \end{bmatrix}

然后我们可以计算仅包含旋转的 Σ\mathbf{\Sigma}

Σ=Σ^L\mathbf{\Sigma} = \mathbf{\hat{\Sigma}}\mathbf{L}

若矩阵中不存在反射，则 L=I\mathbf{L}=\mathbf{I}，所以相当于不做任何操作
但是在计算 U\mathbf{U} 和 V\mathbf{V} 时我们必须小心一点，我们必须找到反射究竟潜藏在 U^\mathbf{\hat{U}} 和 V^\mathbf{\hat{V}} 的哪一个中，并只对那一个矩阵施加反射矩阵 L\mathbf{L}

U={U^L若 det⁡U^<0 且 V^>0U^其他情况V={V^L若 det⁡U^>0 且 V^<0V^其他情况\begin{align*} \mathbf{U} = \begin{cases} \mathbf{\hat{U}}\mathbf{L} & \text{若} \ \det\mathbf{\hat{U}}<0 \ \text{且} \ \mathbf{\hat{V}}>0 \\ \mathbf{\hat{U}} & \text{其他情况} \end{cases} \\ \mathbf{V} = \begin{cases} \mathbf{\hat{V}}\mathbf{L} & \text{若} \ \det\mathbf{\hat{U}}>0 \ \text{且} \ \mathbf{\hat{V}}<0 \\ \mathbf{\hat{V}} & \text{其他情况} \end{cases} \end{align*}

如此我们就构造了仅含有旋转的 SVD，我们也就可以得到仅含有旋转的极分解

R=UVTS=VΣVT\begin{align*} \mathbf{R} &= \mathbf{U}\mathbf{V}^T \\ \mathbf{S} &= \mathbf{V}\mathbf{\Sigma}\mathbf{V}^T \end{align*}

代码如下

function [U Sigma V] = svd_rv (F)
  [U Sigma V] = svd(F);
 
  % reflection matrix
  L = eye (3 ,3);
  L(3 ,3) = det(U * V');
 
  % see where to pull the reflection out of
  detU = det(U);
  detV = det(V);
  if (detU < 0 && detV > 0)
    U = U * L;
  elseif (detU > 0 && detV < 0)
    V = V * L;
  end
 
  % push the reflection to the diagonal
  Sigma = Sigma * L;
end
 
function [R S] = polar_decomposition_rv (F)
  [U Sigma V] = svd_rv (F);
  R = U * V';
  S = V * Sigma * V';
end