An orientation of a rigid body can be described by a rotation matrix \( \mathbf{R} \in SO(3) \), but a 3 × 3 matrix has nine entries with six constraints. Euler angles compress this to three parameters (angles) by expressing any rotation as a sequence of three elemental rotations about coordinate axes (either fixed or body-attached).

A generic Euler-angle parameterization can be written as

\[ \mathbf{R}(\alpha,\beta,\gamma) = \mathbf{R}_{a_1}(\alpha)\,\mathbf{R}_{a_2}(\beta)\,\mathbf{R}_{a_3}(\gamma), \]

where each \( \mathbf{R}_{a_i}(\cdot) \) is a basic rotation about axis \( a_i \in \{x,y,z\} \). Common conventions include intrinsic (body-fixed) and extrinsic (space-fixed) sequences, as well as Tait–Bryan angles (rotations about three distinct axes, like yaw–pitch–roll).

From a robotics modeling viewpoint, Euler angles provide:

• Minimal (3D) local coordinates for \( SO(3) \);
• Direct physical interpretation (yaw, pitch, roll);
• Simple constraints for keeping angles in canonical intervals.

However, we will see that any global 3-parameter representation of \( SO(3) \) must have singularities, and Euler angles suffer from non-uniqueness and severe conditioning near these singular configurations.

flowchart TD
  A["Orientation R in SO(3)"] --> B["Choose Euler convention (e.g. ZYX yaw-pitch-roll)"]
  B --> C["Map to angles q = (phi, theta, psi)"]
  C --> D["Compute R(q) = Rz(psi) * Ry(theta) * Rx(phi)"]
  D --> E["Use R(q) in forward kinematics"]
  C --> F["Differentiate q to get q_dot"]
  F --> G["Relate q_dot to angular velocity omega"]
  G --> H["Detect singularities (det T(q) = 0)"]
        

A widely used convention in robotics is the extrinsic ZYX sequence: a yaw rotation about the world z-axis, followed by a pitch about the world y-axis, followed by a roll about the world x-axis. Let \( \psi \) = yaw, \( \theta \) = pitch, \( \phi \) = roll. Then

\[ \mathbf{R}(\phi,\theta,\psi) = \mathbf{R}_z(\psi)\,\mathbf{R}_y(\theta)\,\mathbf{R}_x(\phi). \]

The elemental rotation matrices are

\[ \mathbf{R}_x(\phi) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\phi & -\sin\phi \\ 0 & \sin\phi & \cos\phi \end{bmatrix},\quad \mathbf{R}_y(\theta) = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix}, \]

\[ \mathbf{R}_z(\psi) = \begin{bmatrix} \cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{bmatrix}. \]

Multiplying yields the explicit ZYX yaw–pitch–roll rotation:

\[ \mathbf{R}(\phi,\theta,\psi) = \begin{bmatrix} c_\psi c_\theta & c_\psi s_\theta s_\phi - s_\psi c_\phi & c_\psi s_\theta c_\phi + s_\psi s_\phi \\ s_\psi c_\theta & s_\psi s_\theta s_\phi + c_\psi c_\phi & s_\psi s_\theta c_\phi - c_\psi s_\phi \\ -s_\theta & c_\theta s_\phi & c_\theta c_\phi \end{bmatrix}, \]

where for brevity \( c_\bullet = \cos(\bullet) \) and \( s_\bullet = \sin(\bullet) \).

A common choice of canonical ranges is

\[ -\pi < \phi \le \pi,\quad -\tfrac{\pi}{2} < \theta < \tfrac{\pi}{2},\quad -\pi < \psi \le \pi. \]

Within this domain, every orientation except certain boundary cases has a unique ZYX representation. However, as we will prove, the mapping is not globally injective nor regular on all of \( SO(3) \).

The mapping \( f : \mathbb{R}^3 \to SO(3) \), \( f(\phi,\theta,\psi) = \mathbf{R}(\phi,\theta,\psi) \) is smooth and surjective but not injective. For example, consider the ZYX angles

\[ (\phi,\theta,\psi) = (0,0,0) \quad\text{and}\quad (\phi,\theta,\psi) = (0,0,2\pi). \]

Both correspond to the identity rotation matrix \( \mathbf{R} = \mathbf{I}_3 \). More subtly, when the middle angle reaches a singular value (see next section), different triplets represent the same rotation.

The deeper reason is topological: the manifold \( SO(3) \) is diffeomorphic to the real projective space \( \mathbb{RP}^3 \), which is compact and has nontrivial topology (it is not homeomorphic to \( \mathbb{R}^3 \)). Therefore, there is no global smooth bijection between \( SO(3) \) and \( \mathbb{R}^3 \).

A rigorous statement is:

\[ \text{There is no smooth map } \Phi: SO(3) \to \mathbb{R}^3 \text{ that is a global coordinate chart.} \]

Any 3-parameter representation (including Euler, axis-angle, etc.) must suffer from singularities or discontinuities. Euler angles are convenient, but this topological fact guarantees there will be configurations where their Jacobians lose rank.

Let \( \boldsymbol{q} = [\phi,\theta,\psi]^{\mathsf{T}} \) be the ZYX Euler angles and \( \dot{\boldsymbol{q}} = [\dot{\phi},\dot{\theta},\dot{\psi}]^{\mathsf{T}} \). The (body-frame) angular velocity \( \boldsymbol{\omega} \in \mathbb{R}^3 \) is related to \( \dot{\boldsymbol{q}} \) via a kinematic Jacobian

\[ \boldsymbol{\omega} = \mathbf{T}(\phi,\theta)\,\dot{\boldsymbol{q}}. \]

For the ZYX convention, one convenient form is

\[ \mathbf{T}(\phi,\theta) = \begin{bmatrix} 1 & 0 & -\sin\theta \\ 0 & \cos\phi & \cos\theta\sin\phi \\ 0 & -\sin\phi & \cos\theta\cos\phi \end{bmatrix}. \]

Proof sketch. Consider an infinitesimal change of angles \( \delta\boldsymbol{q} \). The corresponding incremental rotation matrix is

\[ \delta\mathbf{R} \approx \mathbf{R}_z(\psi)\,\mathbf{R}_y(\theta)\,\left( \tfrac{\partial \mathbf{R}_x(\phi)}{\partial \phi}\,\delta\phi \right) + \mathbf{R}_z(\psi)\,\left( \tfrac{\partial \mathbf{R}_y(\theta)}{\partial \theta}\,\delta\theta \right)\mathbf{R}_x(\phi) + \left( \tfrac{\partial \mathbf{R}_z(\psi)}{\partial \psi}\,\delta\psi \right)\mathbf{R}_y(\theta)\,\mathbf{R}_x(\phi). \]

Using the identity \( \widehat{\boldsymbol{\omega}} = \dot{\mathbf{R}}\mathbf{R}^{\mathsf{T}} \) (where \( \widehat{\boldsymbol{\omega}} \) is the skew-symmetric matrix of \( \boldsymbol{\omega} \)) and matching the coefficients of \( \dot{\phi},\dot{\theta},\dot{\psi} \), one obtains the above \( \mathbf{T}(\phi,\theta) \).

The determinant of the Jacobian is

\[ \det\mathbf{T}(\phi,\theta) = \cos\theta. \]

Thus, whenever \( \cos\theta = 0 \), i.e. \( \theta = \pm \tfrac{\pi}{2} \), the matrix \( \mathbf{T} \) loses rank and cannot be inverted. In such configurations:

• The mapping \( \dot{\boldsymbol{q}} \mapsto \boldsymbol{\omega} \) is still well-defined, but not one-to-one;
• Small changes in orientation may require large changes in Euler angles;
• Control laws expressed in Euler angles may become numerically unstable.

This is the classical gimbal lock phenomenon: at \( \theta = \pm \tfrac{\pi}{2} \), two of the rotation axes line up and the representation effectively loses one degree of freedom in its velocity parameterization.

At a singular configuration (for ZYX: \( \theta = \pm \tfrac{\pi}{2} \)), distinct Euler angle triplets correspond to the same \( \mathbf{R} \). For instance, consider

\[ \theta = \tfrac{\pi}{2},\quad \phi = \alpha,\quad \psi = \beta \]

and

\[ \theta' = \tfrac{\pi}{2},\quad \phi' = \alpha + \delta,\quad \psi' = \beta - \delta. \]

One can verify directly that \( \mathbf{R}(\phi,\theta,\psi) = \mathbf{R}(\phi',\theta',\psi') \) for all \( \delta \). Intuitively, rotating about one axis can be compensated by an opposite rotation about another axis once the gimbals are aligned.

Numerically, near singularities:

• The inverse mapping \( \mathbf{R} \mapsto (\phi,\theta,\psi) \) becomes ill-conditioned;
• Finite-precision arithmetic amplifies noise in angle recovery;
• Integration of Euler angle rates may drift into discontinuous branches unless carefully wrapped back into canonical ranges.

A typical software strategy is to:

1. Detect when \( |\cos\theta| \) falls below a threshold;
2. Switch to an alternative convention, or temporarily represent orientation by a rotation matrix or (later in this chapter) by quaternions;
3. Keep Euler angles only for user interaction and visualization.

flowchart TD
  S["Given Euler angles \n(phi, theta, psi)"] --> C1{"Is abs(cos(theta)) < eps?"}
  C1 -->|no| N["Safe region: \nuse Jacobian T(phi, theta)"]
  C1 -->|yes| G["Near gimbal lock: \nill-conditioned"]
  G --> A1["Option: switch axis \nconvention locally"]
  G --> A2["Option: use rotation matrix / \nquaternion internally"]
        

In this section, we implement ZYX Euler to rotation matrix and the inverse mapping in several languages commonly used in robotics. We also illustrate how to compute angular velocity using \( \mathbf{T}(\phi,\theta) \).

6.1 Python (NumPy and SciPy)

import numpy as np

def euler_zyx_to_R(phi, theta, psi):
    """
    ZYX (yaw-pitch-roll) Euler angles to rotation matrix.
    Angles in radians.
    """
    cphi, sphi = np.cos(phi), np.sin(phi)
    cth,  sth  = np.cos(theta), np.sin(theta)
    cpsi, spsi = np.cos(psi), np.sin(psi)

    R = np.array([
        [cpsi * cth,  cpsi * sth * sphi - spsi * cphi,  cpsi * sth * cphi + spsi * sphi],
        [spsi * cth,  spsi * sth * sphi + cpsi * cphi,  spsi * sth * cphi - cpsi * sphi],
        [-sth,        cth * sphi,                      cth * cphi]
    ])
    return R

def R_to_euler_zyx(R):
    """
    Inverse mapping R -> (phi, theta, psi) for ZYX convention.
    Assumes R is a valid rotation matrix.
    """
    # Protect against numerical issues by clipping
    sth = -R[2, 0]
    sth = np.clip(sth, -1.0, 1.0)
    theta = np.arcsin(sth)

    cth = np.cos(theta)
    if abs(cth) > 1e-6:
        phi = np.arctan2(R[2, 1], R[2, 2])
        psi = np.arctan2(R[1, 0], R[0, 0])
    else:
        # Gimbal lock: theta ~ +/- pi/2
        # psi and phi not uniquely defined; choose one convention
        phi = 0.0
        psi = np.arctan2(-R[0, 1], R[1, 1])

    return phi, theta, psi

def euler_zyx_to_omega(phi, theta, psi, phi_dot, theta_dot, psi_dot):
    """
    Compute body angular velocity omega given Euler angles and their rates.
    """
    cphi, sphi = np.cos(phi), np.sin(phi)
    cth,  sth  = np.cos(theta), np.sin(theta)

    T = np.array([
        [1.0, 0.0, -sth],
        [0.0, cphi, cth * sphi],
        [0.0, -sphi, cth * cphi]
    ])
    qdot = np.array([phi_dot, theta_dot, psi_dot])
    return T @ qdot

# Example usage
if __name__ == "__main__":
    phi, theta, psi = 0.2, -0.3, 0.4
    R = euler_zyx_to_R(phi, theta, psi)
    phi_rec, theta_rec, psi_rec = R_to_euler_zyx(R)
    print("R =\n", R)
    print("Recovered angles:", phi_rec, theta_rec, psi_rec)
      

Robotics-oriented Python libraries such as scipy.spatial.transform.Rotation, pytransform3d, or the modern_robotics toolbox provide tested implementations that internally handle numerical edge cases.

6.2 C++ with Eigen

#include <Eigen/Dense>
#include <cmath>

Eigen::Matrix3d eulerZYXToR(double phi, double theta, double psi) {
    double cphi = std::cos(phi), sphi = std::sin(phi);
    double cth  = std::cos(theta), sth = std::sin(theta);
    double cpsi = std::cos(psi), spsi = std::sin(psi);

    Eigen::Matrix3d R;
    R <<
        cpsi * cth,  cpsi * sth * sphi - spsi * cphi,  cpsi * sth * cphi + spsi * sphi,
        spsi * cth,  spsi * sth * sphi + cpsi * cphi,  spsi * sth * cphi - cpsi * sphi,
        -sth,        cth * sphi,                      cth * cphi;
    return R;
}

Eigen::Vector3d eulerZYXToOmega(double phi, double theta,
                                double phi_dot, double theta_dot, double psi_dot) {
    double cphi = std::cos(phi), sphi = std::sin(phi);
    double cth  = std::cos(theta), sth = std::sin(theta);

    Eigen::Matrix3d T;
    T <<
        1.0, 0.0, -sth,
        0.0, cphi, cth * sphi,
        0.0, -sphi, cth * cphi;

    Eigen::Vector3d qdot(phi_dot, theta_dot, psi_dot);
    return T * qdot;
}
      

In larger robotics codebases, libraries such as KDL, Sophus, kindr, or ROS TF use similar parameterizations internally, often exposing Euler-angle helpers only for configuration interfaces and logging.

6.3 Java Implementation

public final class EulerZYX {

    public static double[][] eulerToR(double phi, double theta, double psi) {
        double cphi = Math.cos(phi), sphi = Math.sin(phi);
        double cth  = Math.cos(theta), sth = Math.sin(theta);
        double cpsi = Math.cos(psi), spsi = Math.sin(psi);

        double[][] R = new double[3][3];
        R[0][0] = cpsi * cth;
        R[0][1] = cpsi * sth * sphi - spsi * cphi;
        R[0][2] = cpsi * sth * cphi + spsi * sphi;

        R[1][0] = spsi * cth;
        R[1][1] = spsi * sth * sphi + cpsi * cphi;
        R[1][2] = spsi * sth * cphi - cpsi * sphi;

        R[2][0] = -sth;
        R[2][1] = cth * sphi;
        R[2][2] = cth * cphi;

        return R;
    }

    public static double[] eulerToOmega(double phi, double theta,
                                        double phiDot, double thetaDot, double psiDot) {
        double cphi = Math.cos(phi), sphi = Math.sin(phi);
        double cth  = Math.cos(theta), sth = Math.sin(theta);

        double[][] T = new double[3][3];
        T[0][0] = 1.0; T[0][1] = 0.0;   T[0][2] = -sth;
        T[1][0] = 0.0; T[1][1] = cphi;  T[1][2] = cth * sphi;
        T[2][0] = 0.0; T[2][1] = -sphi; T[2][2] = cth * cphi;

        double[] qdot = new double[] {phiDot, thetaDot, psiDot};
        double[] omega = new double[3];
        for (int i = 0; i < 3; ++i) {
            omega[i] = 0.0;
            for (int j = 0; j < 3; ++j) {
                omega[i] += T[i][j] * qdot[j];
            }
        }
        return omega;
    }
}
      

In JVM-based robotics stacks, this logic can be wrapped in higher-level classes representing rigid transforms, often alongside matrix and quaternion utilities.

6.4 MATLAB / Simulink

function R = eulerZYXToR(phi, theta, psi)
% ZYX Euler angles to rotation matrix

cphi = cos(phi); sphi = sin(phi);
cth  = cos(theta); sth = sin(theta);
cpsi = cos(psi); spsi = sin(psi);

R = [ cpsi * cth,  cpsi * sth * sphi - spsi * cphi,  cpsi * sth * cphi + spsi * sphi;
      spsi * cth,  spsi * sth * sphi + cpsi * cphi,  spsi * sth * cphi - cpsi * sphi;
      -sth,        cth * sphi,                      cth * cphi ];
end

function omega = eulerZYXToOmega(phi, theta, phi_dot, theta_dot, psi_dot)
% Angular velocity for ZYX Euler angles

cphi = cos(phi); sphi = sin(phi);
cth  = cos(theta); sth = sin(theta);

T = [1,    0,      -sth;
     0,    cphi,   cth * sphi;
     0,   -sphi,   cth * cphi];

omega = T * [phi_dot; theta_dot; psi_dot];
end
      

In Simulink, these functions can be called from a MATLAB Function block to convert state variables (Euler angles) to rotation matrices or angular velocities for kinematic models. The Robotics System Toolbox also provides built-in functions such as eul2rotm and rotm2eul with configurable conventions.

6.5 Wolfram Mathematica

(* ZYX Euler to rotation matrix *)
eulerZYXToR[{phi_, theta_, psi_}] := Module[
  {Rx, Ry, Rz},
  Rx = { {1, 0, 0},
        {0, Cos[phi], -Sin[phi]},
        {0, Sin[phi],  Cos[phi]} };
  Ry = { {Cos[theta], 0, Sin[theta]},
        {0, 1, 0},
        {-Sin[theta], 0, Cos[theta]} };
  Rz = { {Cos[psi], -Sin[psi], 0},
        {Sin[psi],  Cos[psi], 0},
        {0,         0,        1} };
  Rz . Ry . Rx
];

(* Using built-in RotationMatrix: convention ZYX yaw-pitch-roll *)
eulerZYXToR2[{phi_, theta_, psi_}] :=
  RotationMatrix[{psi, theta, phi}, {3, 2, 1}];

(* Inverse: from rotation matrix to Euler angles (ZYX) *)
RToEulerZYX[R_?MatrixQ] := Module[
  {theta, phi, psi, sth, cth},
  sth = -R[[3, 1]];
  theta = ArcSin[Clip[sth, {-1, 1}]];
  cth = Cos[theta];
  If[Abs[cth] > 10^-6,
    phi = ArcTan[R[[3, 2]], R[[3, 3]]];
    psi = ArcTan[R[[2, 1]], R[[1, 1]]],
    (* gimbal lock case *)
    phi = 0.0;
    psi = ArcTan[-R[[1, 2]], R[[2, 2]]]
  ];
  {phi, theta, psi}
];
      

Symbolic manipulation of the rotation matrices in Mathematica is convenient for deriving Jacobians, verifying orthogonality, and exploring singularities.

In this lesson we introduced Euler angles as a 3-parameter representation of orientation built from sequences of basic rotations. Using the ZYX yaw–pitch–roll convention, we derived explicit formulas for the rotation matrix and the relationship between Euler angle rates and angular velocity via the Jacobian \( \mathbf{T}(\phi,\theta) \), whose determinant \( \cos\theta \) reveals kinematic singularities (gimbal lock) at \( \theta = \pm \tfrac{\pi}{2} \).

We linked these singularities to topological obstructions that prevent any global smooth parametrization of \( SO(3) \) by three unconstrained variables. Finally, we implemented Euler-angle conversions in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, emphasizing numerical care near singular configurations. In subsequent lessons, we will introduce alternative representations (axis–angle and unit quaternions) that mitigate some limitations while introducing their own trade-offs.