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

Faster approach for decomposing a rotation to rotations around arbitrary orthogonal axes

我有一个旋转,我想将它分解成围绕 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;
}

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

确实有一个封闭形式的解决方案。由于轴构成标准正交基 A(每个轴是矩阵的一列),您可以通过将 R 转换为基 [=11] 来分解三个轴上的旋转 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]]

我在 python 中做了同样的事情,发现 @Gilles-PhilippePaillé 的回答非常有用,尽管我不得不调整一些东西,主要是使欧拉角反向。我想我会在这里添加我的 python 版本以供参考,以防它对任何人有帮助。

import numpy as np
from scipy.spatial.transform import Rotation

def normalise(v: np.ndarray) -> np.ndarray:
    """Normalise an array along its final dimension."""
    return v / norm(v, axis=-1, keepdims=True)

# Generate random basis
A = Rotation.from_rotvec(normalise(np.random.random(3)) * np.random.rand() * np.pi).as_matrix()

# Generate random rotation matrix
t0 = np.random.rand() * np.pi
t1 = np.random.rand() * np.pi
t2 = np.random.rand() * np.pi
R = Rotation.from_rotvec(A[:, 0] * t0) * Rotation.from_rotvec(A[:, 1] * t1) * Rotation.from_rotvec(A[:, 2] * t2)
R = R.as_matrix()

# Decompose rotation matrix R into the axes of A
rp = Rotation.from_matrix(A.T @ R @ A)
a3, a2, a1 = rp.as_euler('zyx')
xp = A @ Rotation.from_rotvec(a1 * np.array([1, 0, 0])).as_matrix() @ A.T
yp = A @ Rotation.from_rotvec(a2 * np.array([0, 1, 0])).as_matrix() @ A.T
zp = A @ Rotation.from_rotvec(a3 * np.array([0, 0, 1])).as_matrix() @ A.T

# Test that the generated matrix is equal to 'r' (should give 0)
assert np.allclose(xp @ yp @ zp, R)

# Test that the individual rotations preserve the axes (should give 0)
assert np.allclose(xp @ A[:, 0], A[:, 0])
assert np.allclose(yp @ A[:, 1], A[:, 1])
assert np.allclose(zp @ A[:, 2], A[:, 2])