Gimbal Lock

What in this article?

Gimbal lock is a problem originated from the Euler angle parameterization of the 3D rotation.

For 3-axis mechanical gimbal, gimbal lock is a big issue. The 2 most famous gimbal locks happened in the Apollo 11 and Apollo 13 mission.

Most articles and textbooks only describe the Gimbal lock without concrete math. I’d like to show the gimbal lock mathematically and write a short problem to show the gimbal lock.

I will discuss,
1. What is the Euler angle parameterization?
2. What is the Gimbal lock?
3. Write a simple program to show the Gimbal lock.
4. Show how to find the Gimbal lock for a Euler angle order.

The Rotation Group

Mathematically, a rotation is represented as a linear transformation, in another word, a matrix.

The rotation group is the fancy name for the collection of all rotation matrices. The definite of the rotation group \(SO(3)\) is,

\(\{R : R^T R = I, det(R) = 1\}\)

There are many great books and articles about the SO(3) group. I will not go into detail about it in this article.

There exist many methods to parameterize the rotation matrix. The Euler angle is one of them. The idea is a rotation matrix can be represented as 3 consecutive rotation around x or y or z-axis.

The Axis Angle

Axis angle itself is a method to parameterize the rotation. In fact, any rotation matrix can be represented by axis angle. i.e. how can I rotate the original coordinate axis (defined by \(\text{x-axis} := [1, 0, 0], \text{y-axis} := [0, 1, 0], \text{z-axis} := [0,0,1]\)) around a axis \(w\) for \(\theta\) radius so that the rotated coordinate axis align with the coordinate system by the cols of R.

For example, If we want to rotate the original coordinate system around the x axis \(x = [1, 0, 0]\) for \(\theta\) radius, the rotation matrix is defined as,

\(R_x = \begin{bmatrix}
1 & 0 & 0\\
0 & \cos \theta & -\sin \theta\\
0 & \sin \theta & \cos \theta \end{bmatrix}\)

Like wise, if we want to rotate the original coordinate system around the y axis \(y = [0, 1, 0]\) for \(\theta\) radius,

\(R_y = \begin{bmatrix}
\cos \theta & 0 & \sin \theta \\
0 & 1 & 0 \\
-\sin \theta & 0 & \cos \theta \end{bmatrix}\)

if we want to rotate the original coordinate system around the z axis \(z = [0, 0, 1]\) for \(\theta\) radius,

\(R_z = \begin{bmatrix}
\cos \theta & – \sin \theta & 0\\
\sin \theta & \cos \theta & 0 \\
0 & 0 & 1 \end{bmatrix}\)

The Euler Angles

The Euler angles are another way to parameterize the rotation matrix. It implicitly defines 3 rotation axis. Then, It breaks a rotation into 3 axis angle rotations. Because these axis are implicitly defined, the parameters are the magnitude of these 3 rotations.

The most famous Euler angles are the (roll, pitch, yaw) where a rotation is decomposed by first rotating around the z-axis, then around the y-axis, then the x-axis. We call it the ZYX Euler angles. It’s defined as,

\(R_{zyx}(a, b, c) = R_x(a) R_y(b) R_z (c) \; \; \; \; (1)\)

Where (a = roll, b = pitch, c = yaw), and \(R_x(.), R_y(.), R_z(.)\) are defined in the previous section.

(1) is a function which maps \(R^3 \to SO(3)\). But, It is not a Bijective function. For “some” \(x \in R^3\), they are mapped to the same \(y \in SO(3)\).

It’s the Gimbal Lock.

The Gimbal Lock

For the ZYX Eular angles, when \(b = \pm \pi / 2\), there exist infinite number of \((a,b,c)\) which maps to the same \(SO(3)\).

Let’s prove it. Start with the assumption \(b = \pi / 2\),

\(R_{zyx}(a,b,c) = R_x(a) R_y(\pi / 2) R_z(b)\)

\(= \begin{bmatrix}
1 & 0 & 0\\
0 & \cos a & -\sin a\\
0 & \sin a & \cos a \end{bmatrix}

\begin{bmatrix}
0 & 0 & 1\\
0 & 1 & 0\\
-1 & 0 & 0 \end{bmatrix}

\begin{bmatrix}
\cos c & – \sin c & 0\\
\sin c & \cos c & 0 \\
0 & 0 & 1 \end{bmatrix} \\
=
\begin{bmatrix}
1 & 0 & 0\\
0 & \cos a & -\sin a\\
0 & \sin a & \cos a \end{bmatrix}

\begin{bmatrix}
0 & 0 & 1\\
\sin c & \cos c & 0 \\
-\cos c & \sin c & 0 \end{bmatrix} \\
=
\begin{bmatrix}
0 & 0 & 1\\
\cos a \sin c + \sin a \cos c & \cos a \cos c – \sin a \sin c & 0\\
\sin a \sin c – \cos a \cos c & \cos a \sin c + \cos a \sin c & 0 \end{bmatrix}
\)

Recall two useful trigonometric identity. \(\cos a \sin c \pm \sin a \cos c = sin(a \pm c)\) and \(\cos a \cos c \mp \sin a \sin c = cos(a \pm c)\).

\(=
\begin{bmatrix}
0 & 0 & 1\\
\sin (a + c) & \cos (a + c) & 0 \\
-\cos (a+c) & \sin (a + c) & 0 \end{bmatrix}
\\ =
\begin{bmatrix}
0 & 0 & 1\\
\sin d & \cos d & 0 \\
-\cos d & \sin d & 0 \end{bmatrix}
\)

where \(d = a + c\).

This shows when \(a + c = d\), \(R_{zyx}(a, \pi/2, c)\) maps to the same rotation matrix. i.e. when \(b = \pi /2\), each rotation matrix has infinite number of corresponding \((a,c)\).

For a great geometric explanation about gimbal lock please watch this video.

For a mechanical gimbal, the roll, pitch, yaw corresponds to actually physical orientation of the gimbal axis. What will happen if the mechanical gimbal reaches the gimbal lock? Check this video.

Plot the Gimbal Lock

Since we computed the rotation matrix when the gimbal lock happens, I’d like to see how the rotation looks like and compute its rotation axis.

The rotation matrix for gimbal lock is,

\(R(d) = \begin{bmatrix}
0 & 0 & 1\\
\sin d & \cos d & 0 \\
-\cos d & \sin d & 0 \end{bmatrix}\)

To compute the frame after this rotation, we can simple compute,

\(\text{x-new} = R(d) [1, 0, 0]^T\)
\(\text{y-new} = R(d) [0, 1, 0]^T\)
\(\text{z-new} = R(d) [0, 0, 1]^T\)

To compute the rotational axis corresponding to the rotation, we can compute its matrix log and do an unskew. But it is easier to solve for the axis directly. Assuming \(n\) is the axis, by definition of the rotational axis, a point on the axis doesn’t change after applying the rotation. So we have,

\(R(d) n = n\)

\((R(d) – I) n = 0\)

It’s the famous \(Ax = 0\) problem which is ubiquitous in computer vision. It can be solved by SVD. Read these slides for more details.

Put everything together, we can plot the rotated frame and the rotational axis.

Interestingly, \(R(d)\) only has 1 degree of freedom. I expected the rotational axis corresponding to the \(R(d)\) to be “regular”. But the axis is totally irregular.

If you have insight about this problem, please let me know.

Code is here for reference.

import numpy as np
from math import cos, sin
from scipy.linalg import logm, expm

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d


class Arrow3D(FancyArrowPatch):
    ...

def rotate_lock(a):
    r = np.array([
        [0, 0, 1],
        [sin(a), cos(a), 0],
        [-cos(a), sin(a), 0],
    ])

    return r


def solve_axis_svd(R):
    A = R - np.identity(3)
    u, s, v = np.linalg.svd(A)
    n = v[2, :].T

    direction = np.array([1, 0, 0.])
    if n @ direction < 0:
        n = - n

    return n


def solve_axis_log(R):
    w_skew = logm(R)

    w1 = w_skew[2, 1]
    w2 = w_skew[0, 2]
    w3 = w_skew[1, 0]
    return np.array([w1, w2, w3])


def try_plot(R, fig):
    fig.clf()

    ax1 = fig.add_subplot(111, projection='3d')
   
    ...
    
    # Here we create the arrows:
    arrow_prop_dict = dict(
        mutation_scale=20, arrowstyle='->', shrinkA=0, shrinkB=0)

    x1, y1, z1 = R[:, 0]
    a = Arrow3D([0, x1], [0, y1], [0, z1], **arrow_prop_dict, color='r')
    ax1.add_artist(a)

    x2, y2, z2 = R[:, 1]
    a = Arrow3D([0, x2], [0, y2], [0, z2], **arrow_prop_dict, color='b')
    ax1.add_artist(a)

    x3, y3, z3 = R[:, 2]
    a = Arrow3D([0, x3], [0, y3], [0, z3], **arrow_prop_dict, color='g')
    ax1.add_artist(a)

    # Give them a name:
    ax1.text(0.0, 0.0, -0.1, r'$0$')
    ax1.text(x1, y1, z1, r'$x-new$')
    ax1.text(x2, y2, z2, r'$y-new$')
    ax1.text(x3, y3, z3, r'$z-new$')

    # ref
    a = Arrow3D([0, 1], [0, 0], [0, 0], **
                arrow_prop_dict, color='r', alpha=0.3)
    ax1.add_artist(a)

    a = Arrow3D([0, 0], [0, 1], [0, 0], **
                arrow_prop_dict, color='b', alpha=0.3)
    ax1.add_artist(a)

    a = Arrow3D([0, 0], [0, 0], [0, 1], **
                arrow_prop_dict, color='g', alpha=0.3)
    ax1.add_artist(a)

    # Give them a name:
    ax1.text(1.1, 0, 0, r'$x-ref$')
    ax1.text(0, 1.1, 0, r'$y-ref$')
    ax1.text(0, 0, 1.1, r'$z-ref$')

    # rotation_axis = solve_axis_log(R)
    rotation_axis = solve_axis_svd(R)

    w1, w2, w3 = rotation_axis
    a = Arrow3D([0, w1], [0, w2], [0, w3], **
                arrow_prop_dict, color='m', alpha=0.5)
    ax1.text(w1, w2, w3, r'$w$')
    ax1.add_artist(a)
    # plt.show()


def main():
    fig = plt.figure()
    for i, a in enumerate(np.linspace(0, 2 * np.pi, 60)):
        R = rotate_x(a)
        R = rotate_lock(a)
        try_plot(R, fig)
        plt.pause(0.1)
        plt.savefig('res/{0:03d}.png'.format(i))


if __name__ == "__main__":
    main()

How to compute the gimbal lock for Euler angles?

We know for ZYX Euler angles, when pitch is \(\pm \pi / 2\), we will have the gimbal lock.

An interesting question that came to my mind is: For a Euler angles sequence, how to compute the gimbal lock pose for it?

Before that, we need some more background about the Euler angles. Euler angles can be divided into 2 categories.

1. The Proper Euler angles (z–x–z, x–y–x, y–z–y, z–y–z, x–z–x, y–x–y)
2. The Tait–Bryan angles (x–y–z, y–z–x, z–x–y, x–z–y, z–y–x, y–x–z)

For the Proper Euler angles, it’s easy to see when the middle orientation is 0, we have a gimbal lock.

Consider an example. for the ZXZ Euler angles, when the rotation around X-axis is 0, the first rotation around Z and the last rotation around Z collides.

For the Tait-Bryan angles, singularity happens when, in the global frame, the second rotation makes the last rotational axis aligns with the first rotational axis.

Consider the ZYX Euler angle, since the first rotation \(R_z\) doesn’t change the z-axis in the second rotational frame, to compute the second rotation \(R_y\) which makes the last rotation axis (the x-axis of the current frame) align with its z-axis, we just need to solve,

\(
\begin{bmatrix}
\cos \theta & 0 & \sin \theta \\
0 & 1 & 0 \\
-\sin \theta & 0 & \cos \theta \end{bmatrix}
\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}
=
\begin{bmatrix} 0 \\ 0 \\ \mp 1 \end{bmatrix}\)

It can be simplified to,

\(\begin{bmatrix} \cos \theta = 0 \\ \sin \theta = \pm 1 \end{bmatrix}\)

Solve it we get the gimbal lock for ZYX which is \(R_y(\pm \pi / 2)\).

Let’s do another case.

Consider the YZX Euler angle, following the same procedure (how the second rotation \(R_z\) makes its x-axis align with its y-axis), we got this equation,

\(R_z = \begin{bmatrix}
\cos \theta & – \sin \theta & 0\\
\sin \theta & \cos \theta & 0 \\
0 & 0 & 1 \end{bmatrix}
\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}
=
\begin{bmatrix} 0 \\ \pm 1 \\ 0 \end{bmatrix}\)

It can be simplified to,

\(\begin{bmatrix} \cos \theta = 0 \\ \sin \theta = \pm 1 \end{bmatrix}\)

Solve it we get the gimbal lock for YZX which is \(R_z(\pm \pi / 2)\).

We can show for the Tait–Bryan angles, the gimbal lock happens when the middle rotation is \(\pm \pi / 2\).

(You can also “feel” this result by rotating your finger according to the Euler angle.)

Better Parameterization

Euler angle is a great rotational parameterization because it’s human-readable. However, for the most computational problem, it’s not a good rotational parameterization.

Consider an optimization problem with rotation as its parameter.

\(\text{minimize} \; \; \text{cost}(R)\)

If we use the Euler angle ZYX parameterization, at the Gimbal lock, the cost lost 1 degree of freedom. As a result, the optimization is “unconstrained” by the Euler angle parameterization and the optimization will be slow. We can also argue that when the rotation close to the Gimbal lock, the gradient w.r.t yaw and roll become very dependent. And the optimization is less effective.

Computationally, the mathematical tool to handing this case it the axis-angle parameterization. Please read this article for a simple introduction to axis-angle parameterization and optimization using Lie-group.