将旋转分解为围绕任意正交轴的旋转的更快方法

Question

我有一个旋转，我想将它分解成一系列围绕 3 个正交任意轴的旋转。这有点像欧拉分解的推广，其中旋转不围绕 X、Y 和 Z 轴进行

我试图找到一个封闭形式的解决方案，但没有成功，所以我生成了一个数值解决方案，该解决方案基于最小化我想要的旋转与代表 3 个角度的 3 个轴旋转的 3 个四元数的乘积之间的差异未知数。 'SimplexMinimize' 只是代码的抽象，用于找到最小化错误的 3 个角度。

double GSUtil::ThreeAxisDecomposition(const Quaternion &target, const Vector &ax1, const Vector &ax2, const Vector &ax3, double *ang1, double *ang2, double *ang3)
{
    DataContainer data = {target, ax1, ax2, ax3};
    VaraiablesContainer variables = {ang1, ang2, ang3};
    error = SimplexMinimize(ThreeAxisDecompositionError, data, variables);
}

double GSUtil::ThreeAxisDecompositionError(const Quaternion &target, const Vector &ax1, const Vector &ax2, const Vector &ax3, double ang1, double ang2, double ang3)
{
    Quaternion product = MakeQFromAxisAngle(ax3, ang3) * MakeQFromAxisAngle(ax2, ang2) * MakeQFromAxisAngle(ax1, ang1);
    // now we need a distance metric between product and target. I could just calculate the angle between them:
    // theta = acos(2?q1,q2?^2-1) where ?q1,q2? is the inner product (n1n2 + x1x2+ y1y2 + z1z2)
    // but there are other quantities that will do a similar job in less time
    // 1-(q1,q2)^2 should be faster to calculate and is 0 when they are identical and 1 when they are 180 degrees apart
    double innerProduct = target.n * product.n + target.v.x * product.v.x + target.v.x * product.v.x + target.v.x * product.v.x;
    double error = 1 - innerProduct * innerProduct;
    return error;
}

它有效（我认为）但显然它很慢。我的感觉是应该有一个封闭形式的解决方案。至少该函数应该有一个渐变，这样我就可以使用更快的优化器。

Answer 1

确实有一个封闭形式的解决方案。由于坐标轴形成正交基A（每个坐标轴是矩阵的一列），您可以通过将R 转换为基A 来分解R 在三个坐标轴上的旋转，然后进行欧拉角分解在三个主轴上：

R = A*R'*A^t = A*X*Y*Z*A^t = (A*X*A^t)*(A*Y*A^t)*(A*Z*A^t)

这转化为以下算法：

• 计算R' = A^t*R*A
• 围绕主轴将R'分解为欧拉角以获得矩阵X, Y, Z
• 计算围绕给定轴的三个旋转：
  • X' = A*X*A^t
  • Y' = A*Y*A^t
  • Z' = A*Y*A^t

作为参考，这是我用来测试答案的 Mathematica 代码

(*Generate random axes and a rotation matrix for testing purposes*)
a = RotationMatrix[RandomReal[{0, \[Pi]}], 
   Normalize[RandomReal[{-1, 1}, 3]]];
t1 = RandomReal[{0, \[Pi]}];
t2 = RandomReal[{0, \[Pi]}];
t3 = RandomReal[{0, \[Pi]}];
r = RotationMatrix[t1, a[[All, 1]]].
    RotationMatrix[t2, a[[All, 2]]].
    RotationMatrix[t2, a[[All, 3]]];

(*Decompose rotation matrix 'r' into the axes of 'a'*)
rp = Transpose[a].r.a;
{a1, a2, a3} = EulerAngles[rp, {1, 2, 3}];
xp = a.RotationMatrix[a1, {1, 0, 0}].Transpose[a];
yp = a.RotationMatrix[a2, {0, 1, 0}].Transpose[a];
zp = a.RotationMatrix[a3, {0, 0, 1}].Transpose[a];

(*Test that the generated matrix is equal to 'r' (should give 0)*)
xp.yp.zp - r // MatrixForm

(*Test that the individual rotations preserve the axes (should give 0)*)
xp.a[[All, 1]] - a[[All, 1]]
yp.a[[All, 2]] - a[[All, 2]]
zp.a[[All, 3]] - a[[All, 3]]