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

Question

I have a rotation and I want to decompose it into a series of rotations around 3 orthogonal arbitrary axes. It's a bit like a generalisation of Euler decomposition where the rotations are not around the X, Y and Z axes

I've tried to find a closed form solution but not been successful so I have produced a numerical solution based on minimising the difference between the rotation I want and the product of 3 quaternions representing the 3 axes roations with the 3 angles being the unknowns. 'SimplexMinimize' is just an abstraction of the code to find the 3 angles that minimises the error.

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;
}

It works (I think) but obviously it is quite slow. My feeling is there ought to be a closed form solution. At the very least there ought to be a gradient to the function so I can use a faster optimiser.

Answer 1

There is indeed a closed form solution. Since the axes form an orthonormal basis A (each axe is a column of the matrix), you can decompose a rotation R on the three axes by transforming R into the basis A and then do Euler Angle decomposition on the three main axes:

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

This translates into the following algorithm:

• Compute R' = A^t*R*A
• Decompose R' into Euler Angles around main axes to obtain matrices X, Y, Z
• Compute the three rotations around the given axes:
  • X' = A*X*A^t
  • Y' = A*Y*A^t
  • Z' = A*Y*A^t

As a reference, here's the Mathematica code I used to test my answer

(*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]]

Answer 2

I was doing the same thing in python and found @Gilles-PhilippePaillé 's answer really helpful although I had to tweak a couple of things, mostly getting the euler angles out in reverse. Thought I would add my python version here for reference anyway in case it helps anyone.

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])