2

I have a $4\times 4$ matrix $M$ which represents a general 4-dimensional rotation. $$ M = \pmatrix{a_{11} &a_{12} &a_{13} &a_{14}\\a_{21} &a_{22} &a_{23} &a_{24}\\a_{31} &a_{32} &a_{33} &a_{34}\\a_{41} &a_{42} &a_{43} &a_{44}} $$ There are also six 4D basis rotation matrices $M_{XY}$, $M_{YZ}$, $M_{ZX}$, $M_{XW}$, $M_{YW}$, $M_{ZW}$.

How to decompose 4D rotation $M$ into basis rotations (i.e., to find corresponding angles $T_{XY}$, $T_{YZ}$, $T_{ZX}$, $T_{XW}$, $T_{YW}$, $T_{ZW}$)?

Wikipedia says that any 4D rotation given by its matrix can be decomposed into two isoclinic rotations. However, I do not know if it is possible to decompose them further into basis rotations. Also, there may be another (simpler) way.

3 Answers3

1

You should be able to build your basis matrices by doing something similar to the Gram-Schmidt process: by multiplying by suitably-chosen rotation matrices you can set each of the off-diagonal coefficients in turn equal to 0. For instance, consider applying a $M_{ZW}$ rotation by some angle $\theta$ to the matrix: $$M' = M_{ZW}(\theta)\ M = \pmatrix{1&0&0&0\\0&1&0&0\\0&0&\cos\theta&-\sin\theta\\0&0&\sin\theta&\cos\theta}\pmatrix{a_{11}&a_{12}&a_{13}&a_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\a_{31}&a_{32}&a_{33}&a_{34}\\a_{41}&a_{42}&a_{43}&a_{44}} =\pmatrix{a_{11}&a_{12}&a_{13}&a_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\*&*&a_{33}\cos\theta-a_{43}\sin\theta&a_{34}\cos\theta-a_{44}\sin\theta\\*&*&a_{33}\sin\theta+a_{43}\cos\theta&a_{34}\sin\theta+a_{44}\cos\theta}$$ (Here the *s are 'don't-cares' - they have a value, but it's essentially unimportant.) Now, by choosing $\theta$ correctly, you can make $M'_{34} = a_{34}\cos\theta-a_{44}\sin\theta$ equal to zero without affecting any other upper-triangular elements. Once you've got that matrix, you can continue in a similar fashion by finding a suitable $M_{YW}(\phi)$ to zero the $M_{24}$ element, then an $M_{XW}$ to zero the $M_{14}$ element, etc.; repeating this procedure for each combination of unequal coordinates (in some canonically-defined order) will eventually lead you to an orthogonal (by hypothesis) lower-triangular matrix, which must therefore be the identity.

Steven Stadnicki
  • 53,147
0

If what you are looking for is a way to decompose a 4x4 rotation matrix into 4D Euler angles, here is a C++ function that performs this operation (license: The Unlicense / Public Domain). The C++ code:

Euler4D Euler4D::from_basis(const Basis4D &p_basis, const bool p_use_special_cases) {
    // This process produces a LOT of floating-point error.
    // It's also likely that this could be optimized a lot.
    // With 32-bit floats, this has a margin of error of 0.1%.

Basis4D basis = p_basis.orthonormalized();

// Stage 0: Determine if this is a special case. This is optional, for UX.
// This ensures that single rotations in the inspector don't transform into
// other rotations. Godot does the same thing for 3D euler angles.
if (p_use_special_cases) {
    if (Math::is_zero_approx(basis.y.z) && Math::is_zero_approx(basis.z.y) && Math::is_zero_approx(basis.z.x) && Math::is_zero_approx(basis.x.z) && Math::is_zero_approx(basis.x.w) && Math::is_zero_approx(basis.w.x) && Math::is_zero_approx(basis.w.y) && Math::is_zero_approx(basis.y.w)) {
        // Special case: No YZ, ZX, XW, or WY rotations, so the only rotations may be in XY or ZW.
        return Euler4D(0.0, 0.0, Math::atan2(basis.x.y, basis.x.x), 0.0, 0.0, Math::atan2(basis.z.w, basis.z.z));
    }
    if (Math::is_zero_approx(basis.z.x) && Math::is_zero_approx(basis.x.z) && Math::is_zero_approx(basis.x.y) && Math::is_zero_approx(basis.y.x) && Math::is_zero_approx(basis.w.y) && Math::is_zero_approx(basis.y.w) && Math::is_zero_approx(basis.z.w) && Math::is_zero_approx(basis.w.z)) {
        // Special case: No ZX, XY, WY, or ZW rotations, so the only rotations may be in YZ or XW.
        return Euler4D(Math::atan2(basis.y.z, basis.y.y), 0.0, 0.0, Math::atan2(basis.x.w, basis.x.x), 0.0, 0.0);
    }
    if (Math::is_zero_approx(basis.y.z) && Math::is_zero_approx(basis.z.y) && Math::is_zero_approx(basis.x.y) && Math::is_zero_approx(basis.y.x) && Math::is_zero_approx(basis.x.w) && Math::is_zero_approx(basis.w.x) && Math::is_zero_approx(basis.z.w) && Math::is_zero_approx(basis.w.z)) {
        // Special case: No YZ, XY, XW, or ZW rotations, so the only rotations may be in ZX or WY.
        return Euler4D(0.0, Math::atan2(basis.z.x, basis.z.z), 0.0, 0.0, Math::atan2(basis.w.y, basis.w.w), 0.0);
    }
}

// Stage 1: The outer rotations (and yz is special).
real_t yz = -Math::asin(basis.z.y);
real_t zx = Math::atan2(basis.z.x, basis.z.z);
real_t wy = Math::atan2(basis.w.y, basis.y.y);

// Un-rotate by perpendicular ZX and WY, the outermost rotations.
basis = Basis4D::from_zx(-zx) * basis * Basis4D::from_wy(-wy);

// Stage 2: The inner rotations (and xw is special).
real_t xw = -Math::asin(basis.w.x);
real_t zw = Math::atan2(basis.z.w, basis.z.z);
real_t xy = Math::atan2(basis.x.y, basis.y.y);

return Euler4D(yz, zx, xy, xw, wy, zw);



}

Note that this code is designed with the conventions of Y-up, YZ is pitch, and uses ZX and WY planes instead of XZ and YW. Also note that basis components are column-major (basis.x.y is X column Y row).

Note that this is not the same as a generalized 4D geometric algebra rotor, which can describe and interpolate between rotations. This approach heavily depends on the conventions used, cannot be smoothly interpolated, and is prone to gimbal lock.

Aaron Franke
  • 438
-1

You can think of this process in terms of applying successive coordinate-plane rotations, each one on the right of what has gone before, so that ultimately the resulting rotation is the identity.

Then multiplying the identity matrix by the inverses of these coordinate-plane rotations, each one also on the right of what has gone before, you will end up with an equation showing how to express the original rotation as a product of coordinate-plane rotations.

This actually works in any dimension. Since the dimension of SO(n) is n-choose-2 = n(n-1)/2, that should be the maximum number of coordinate-plane rotations needed.

In a given dimension, selecting a fixed ordering for performing this procedure defines a continuous map from the n(n-1)/2 dimensional torus (since each angle may be thought of as a point on circle) to SO(n+1):

Tn(n-1)/2 —> SO(n)

that is in fact surjective (since each rotation can be expressed as a product of coordinate-plane rotations). The circle factors of Tn(n-1)/2 are indexed by pairs $(i, j)$ with $1 \le i \lt j \le n$ (since each one represents one coordinate plane).

The 2-plane rotations used in this procedure are sometimes called Householder rotations.

BLAZE
  • 1
Daz
  • 1