Any rigid-body orientation in 3D can be represented by a rotation matrix \( R \in \mathrm{SO}(3) \), acting on a point \( p \in \mathbb{R}^3 \) as

\[ p_{\text{rot}} = R\,p. \]

The axis–angle representation encodes the same orientation by a unit axis \( \omega \in \mathbb{R}^3 \) and a scalar angle \( \theta \in \mathbb{R} \):

\[ \|\omega\| = 1, \quad R = R(\omega,\theta) \in \mathrm{SO}(3) \quad\text{“rotation by angle } \theta \text{ about axis } \omega\text{”}. \]

The pair \( (\omega,\theta) \) is not unique: for example, \( (\omega,\theta) \) and \( (-\omega,-\theta) \) encode the same rotation, and so do \( (\omega,\theta+2k\pi) \) for any integer \( k \). A common convention is to choose \( 0 \le \theta \le \pi \). When \( \theta = 0 \), the rotation is the identity and the axis is undefined (any unit vector is acceptable).

It is often convenient to combine axis and angle into a rotation vector \( \phi \in \mathbb{R}^3 \):

\[ \phi = \theta\,\omega, \quad \theta = \|\phi\|, \quad \omega = \begin{cases} \phi / \|\phi\| & \text{if } \phi \neq 0,\\[4pt] \text{arbitrary unit vector} & \text{if } \phi = 0. \end{cases} \]

This gives a minimal 3-parameter representation of orientations consistent with the Lie algebra \( \mathfrak{so}(3) \).

flowchart TD
  A["Unit axis omega in R^3"] --> C["Axis-angle pair (omega, theta)"]
  B["Angle theta (radians)"] --> C
  C --> D["Rotation vector phi = theta * omega"]
  D --> E["Rotation matrix R in SO(3) via exp(phi^)^"]
        

Here “phi^” denotes the skew-symmetric matrix corresponding to the vector \( \phi \), introduced formally in the next section.

Let \( \omega = [\omega_1,\omega_2,\omega_3]^\top \in \mathbb{R}^3 \). The hat (or skew-symmetric) operator \( (\cdot)^\wedge : \mathbb{R}^3 \rightarrow \mathfrak{so}(3) \) is defined by

\[ \omega^\wedge \;=\; \begin{bmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{bmatrix}. \]

This matrix satisfies \( \omega^\wedge v = \omega \times v \) for all \( v \in \mathbb{R}^3 \), and every \( X \in \mathfrak{so}(3) \) is of the form \( X = \omega^\wedge \) for a unique vector \( \omega \) (the “vee” operator \( (\cdot)^\vee \) recovers \( \omega \) from \( X \)).

A crucial algebraic identity is obtained by applying \( \omega^\wedge \) twice:

\[ (\omega^\wedge)^2 v = \omega^\wedge(\omega \times v) = \omega(\omega^\top v) - \|\omega\|^2 v, \quad \forall v \in \mathbb{R}^3. \]

Thus, as a matrix,

\[ (\omega^\wedge)^2 = \omega\omega^\top - \|\omega\|^2 I_3, \]

and for a unit axis \( \|\omega\| = 1 \) we have

\[ (\omega^\wedge)^2 = \omega\omega^\top - I_3. \]

2.1 Exponential Map and Rodrigues’ Formula

From Chapter 2, the exponential map from the Lie algebra \( \mathfrak{so}(3) \) to \( \mathrm{SO}(3) \) is given by the matrix exponential. For a unit axis \( \omega \) and angle \( \theta \), define \( K = \omega^\wedge \). Then

\[ R(\omega,\theta) = \exp(\theta K) = I_3 + \theta K + \frac{\theta^2}{2!}K^2 + \frac{\theta^3}{3!}K^3 + \cdots. \]

Using \( K^2 = \omega\omega^\top - I_3 \) and \( K\omega = 0 \), we can show that \( K^3 = -K \) and \( K^4 = -K^2 \), so all higher powers of \( K \) reduce to linear combinations of \( K \) and \( K^2 \). Grouping even and odd powers:

\[ \begin{aligned} R(\omega,\theta) &= I_3 + \left(\theta - \frac{\theta^3}{3!} + \frac{\theta^5}{5!} - \cdots\right)K + \left(\frac{\theta^2}{2!} - \frac{\theta^4}{4!} + \cdots\right)K^2 \\ &= I_3 + \sin\theta \, K + (1-\cos\theta)\,K^2. \end{aligned} \]

This is the Rodrigues rotation formula:

\[ R(\omega,\theta) = I_3 + \sin\theta\,\omega^\wedge + (1-\cos\theta)(\omega^\wedge)^2, \quad \|\omega\| = 1. \]

It gives an efficient way to compute a rotation matrix from an axis–angle pair using only trigonometric functions and matrix multiplications.

Given \( R \in \mathrm{SO}(3) \), we want to recover \( (\omega,\theta) \) such that \( R = \exp(\theta \omega^\wedge) \). The eigenvalues of any rotation matrix are \( 1, e^{\pm i\theta} \), so the trace is

\[ \operatorname{tr}(R) = 1 + 2\cos\theta. \]

This can also be derived directly from Rodrigues’ formula using \( \operatorname{tr}(\omega^\wedge) = 0 \) and \( \operatorname{tr}(\omega\omega^\top) = \|\omega\|^2 = 1 \):

\[ \begin{aligned} \operatorname{tr}(R) &= \operatorname{tr}(I_3) + \sin\theta \operatorname{tr}(\omega^\wedge) + (1-\cos\theta)\operatorname{tr}(\omega\omega^\top - I_3) \\ &= 3 + 0 + (1-\cos\theta)(1-3) \\ &= 3 - 2(1-\cos\theta) = 1 + 2\cos\theta. \end{aligned} \]

Solving for \( \theta \) with the convention \( 0 \le \theta \le \pi \):

\[ \theta = \cos^{-1}\!\left(\frac{\operatorname{tr}(R) - 1}{2}\right), \quad 0 \le \theta \le \pi. \]

3.1 Recovering the Axis

For \( 0 < \theta < \pi \), using \( R = \exp(\theta\omega^\wedge) \) and Rodrigues’ formula, the skew-symmetric part of \( R \) isolates \( \omega^\wedge \):

\[ \begin{aligned} R - R^\top &= \bigl(I_3 + \sin\theta\,\omega^\wedge + (1-\cos\theta)(\omega^\wedge)^2\bigr) \\ &\quad\;\; - \bigl(I_3 - \sin\theta\,\omega^\wedge + (1-\cos\theta)(\omega^\wedge)^2\bigr) \\ &= 2\sin\theta\,\omega^\wedge. \end{aligned} \]

Hence,

\[ \omega^\wedge = \frac{1}{2\sin\theta}(R - R^\top), \quad \omega = \bigl(\omega^\wedge\bigr)^\vee. \]

Writing this out componentwise gives

\[ \omega = \frac{1}{2\sin\theta} \begin{bmatrix} R_{32} - R_{23} \\ R_{13} - R_{31} \\ R_{21} - R_{12} \end{bmatrix}, \quad 0 < \theta < \pi. \]

3.2 Special Cases

• Case \( \theta \approx 0 \). When \( \theta \rightarrow 0 \), \( R \rightarrow I_3 \) and numerical computation of \( \omega \) from \( R-R^\top \) is ill-conditioned. In this limit, any axis is valid and \( \omega \) is typically obtained from application context or from a limiting process.
• Case \( \theta \approx \pi \). The formula using \( R - R^\top \) becomes unreliable because \( \sin\theta \approx 0 \). Instead, one uses the diagonal entries of \( R \) to recover \( \omega \) up to sign, then enforces \( \|\omega\| = 1 \).

Combining these steps defines the matrix logarithm on \( \mathrm{SO}(3) \):

\[ \log(R) = \theta\,\omega^\wedge \in \mathfrak{so}(3), \quad \exp(\log(R)) = R. \]

Define the rotation vector \( \phi \in \mathbb{R}^3 \) by \( \phi = \theta \omega \). Then \( \phi^\wedge = \theta \omega^\wedge \), and the exponential map can be written compactly as

\[ R = \exp(\phi^\wedge), \quad \phi \in \mathbb{R}^3, \; \|\phi\| = \theta. \]

In these coordinates, the origin \( \phi = 0 \) corresponds to the identity rotation. For small \( \|\phi\| \), we have the first-order approximation

\[ R \approx I_3 + \phi^\wedge, \]

which is fundamental for differential kinematics and linearization, and will be used in later chapters when relating orientation errors to twist errors.

The axis–angle pair \( (\omega,\theta) \) and the rotation vector \( \phi \) are equivalent: they contain the same information, but are used in different numerical algorithms depending on which is more convenient.

Implementations must handle round-off errors (e.g., ensuring \( R \) stays in \( \mathrm{SO}(3) \)) and singular cases (\( \theta \approx 0 \) or \( \theta \approx \pi \)). A typical conversion pipeline is summarized below.

flowchart TD
  S["Start"] --> C1{"Input type?"}
  C1 -->|"axis-angle \n(omega, theta)"| F1["Normalize omega; \nbuild skew K"]
  F1 --> F2["Compute R = I + sin(theta) K + \n(1 - cos(theta)) K*K"]
  C1 -->|rotation matrix R| B1["Compute theta from trace(R)"]
  B1 --> B2{"theta near 0 or pi?"}
  B2 -->|yes| B3["Use special-case \nformulas / limits"]
  B2 -->|no| B4["Compute omega from \n(R - R^T)/(2 sin(theta))"]
  F2 --> OUT["Return rotation matrix R"]
  B3 --> OUT2["Return axis-angle \n(omega, theta)"]
  B4 --> OUT2
        

In all numeric code, it is standard to:

• Clamp \( x = (\operatorname{tr}(R) - 1)/2 \) to the interval \( [-1,1] \) before evaluating \( \theta = \cos^{-1}(x) \).
• Use a small threshold \( \varepsilon \) to test \( \theta < \varepsilon \) or \( \pi - \theta < \varepsilon \).
• Optionally re-orthogonalize numerically computed rotations (e.g., by SVD) to project back onto \( \mathrm{SO}(3) \).

Using NumPy, we can implement from-scratch conversions between axis–angle and rotation matrices. This aligns with the Product-of-Exponentials framework used later in the course.

import numpy as np

def skew(omega: np.ndarray) -> np.ndarray:
    """
    Return the skew-symmetric matrix omega^ from a 3-vector omega.
    """
    wx, wy, wz = omega
    return np.array([[0.0, -wz,  wy],
                     [wz,  0.0, -wx],
                     [-wy, wx,  0.0]])

def axis_angle_to_R(omega: np.ndarray, theta: float) -> np.ndarray:
    """
    Rodrigues formula: R = I + sin(theta) K + (1 - cos(theta)) K^2
    where K = skew(omega), ||omega|| = 1.
    """
    omega = np.asarray(omega, dtype=float).reshape(3)
    n = np.linalg.norm(omega)
    if n < 1e-12:
        # Identity rotation
        return np.eye(3)
    omega = omega / n
    K = skew(omega)
    I = np.eye(3)
    return I + np.sin(theta) * K + (1.0 - np.cos(theta)) * (K @ K)

def R_to_axis_angle(R: np.ndarray):
    """
    Inverse map: given R in SO(3), return (omega, theta).
    """
    R = np.asarray(R, dtype=float).reshape(3, 3)
    # Ensure R is close to orthogonal if coming from noisy data
    trace = np.trace(R)
    cos_theta = (trace - 1.0) / 2.0
    cos_theta = np.clip(cos_theta, -1.0, 1.0)
    theta = np.arccos(cos_theta)
    eps = 1e-8

    if theta < eps:
        # Very small rotation: axis arbitrary, choose x-axis
        return np.array([1.0, 0.0, 0.0]), 0.0

    # Use the skew-symmetric part
    omega = np.array([
        R[2, 1] - R[1, 2],
        R[0, 2] - R[2, 0],
        R[1, 0] - R[0, 1]
    ]) / (2.0 * np.sin(theta))

    # Normalize for robustness
    omega = omega / np.linalg.norm(omega)
    return omega, theta

if __name__ == "__main__":
    # Example: rotation of 90 deg about z-axis
    omega = np.array([0.0, 0.0, 1.0])
    theta = np.pi / 2.0
    R = axis_angle_to_R(omega, theta)
    print("R =\n", R)

    omega_rec, theta_rec = R_to_axis_angle(R)
    print("Recovered omega:", omega_rec)
    print("Recovered theta:", theta_rec)
      

In robotics-oriented Python libraries, these operations are encapsulated in high-level classes. For example, with spatialmath and roboticstoolbox-python:

from spatialmath import SO3

omega = [0.0, 0.0, 1.0]
theta = np.pi / 2.0

R = SO3.AngleAxis(theta, omega)        # axis-angle to rotation matrix
omega2, theta2 = R.angvec()            # rotation matrix back to axis-angle
      

These libraries will be helpful when implementing forward kinematics using the PoE formulation in later chapters.

In C++, the Eigen library provides both basic linear algebra and an AngleAxisd class that directly represents axis–angle rotations. Below is a from-scratch implementation compatible with Eigen matrices.

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

using Eigen::Matrix3d;
using Eigen::Vector3d;

Matrix3d skew(const Vector3d& w) {
    Matrix3d K;
    K <<   0.0,   -w.z(),  w.y(),
            w.z(),  0.0,   -w.x(),
           -w.y(),  w.x(),  0.0;
    return K;
}

Matrix3d axisAngleToR(const Vector3d& axis, double theta) {
    Vector3d w = axis;
    double n = w.norm();
    if (n < 1e-12) {
        return Matrix3d::Identity();
    }
    w /= n;
    Matrix3d K = skew(w);
    Matrix3d I = Matrix3d::Identity();
    return I + std::sin(theta) * K + (1.0 - std::cos(theta)) * (K * K);
}

void RToAxisAngle(const Matrix3d& R, Vector3d& axis, double& theta) {
    double trace = R.trace();
    double cos_theta = (trace - 1.0) / 2.0;
    if (cos_theta > 1.0) cos_theta = 1.0;
    if (cos_theta < -1.0) cos_theta = -1.0;
    theta = std::acos(cos_theta);

    const double eps = 1e-8;
    if (theta < eps) {
        axis = Vector3d(1.0, 0.0, 0.0);
        theta = 0.0;
        return;
    }

    axis.x() = (R(2, 1) - R(1, 2)) / (2.0 * std::sin(theta));
    axis.y() = (R(0, 2) - R(2, 0)) / (2.0 * std::sin(theta));
    axis.z() = (R(1, 0) - R(0, 1)) / (2.0 * std::sin(theta));
    axis.normalize();
}
      

Eigen also offers the convenience wrapper:

Eigen::AngleAxisd aa(theta, axis.normalized());
Matrix3d R = aa.toRotationMatrix();
      

In ROS and ROS 2, orientation data is often exchanged as quaternions (geometry_msgs::Quaternion) and transformed using tf2. The tf2::Matrix3x3 and tf2::Quaternion classes internally implement conversions equivalent to the axis–angle formulas derived here.

Java does not have a single standard robotics library, but many projects rely on EJML or Apache Commons Math for linear algebra. Here we show a minimal implementation using plain arrays; the same logic can be wrapped around EJML matrices.

public class AxisAngleSO3 {

    public static double[][] skew(double[] w) {
        double wx = w[0], wy = w[1], wz = w[2];
        return new double[][]{
            { 0.0,  -wz,   wy},
            { wz,    0.0, -wx},
            {-wy,    wx,   0.0}
        };
    }

    public static double[][] axisAngleToR(double[] axis, double theta) {
        double[] w = axis.clone();
        double norm = Math.sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]);
        if (norm < 1e-12) {
            return new double[][]{
                {1.0, 0.0, 0.0},
                {0.0, 1.0, 0.0},
                {0.0, 0.0, 1.0}
            };
        }
        w[0] /= norm; w[1] /= norm; w[2] /= norm;

        double[][] K = skew(w);
        double[][] K2 = matMul(K, K);
        double[][] R = new double[3][3];

        double s = Math.sin(theta);
        double c = Math.cos(theta);

        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                double delta = (i == j) ? 1.0 : 0.0;
                R[i][j] = delta + s * K[i][j] + (1.0 - c) * K2[i][j];
            }
        }
        return R;
    }

    public static double[][] matMul(double[][] A, double[][] B) {
        double[][] C = new double[3][3];
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                double sum = 0.0;
                for (int k = 0; k < 3; ++k) {
                    sum += A[i][k] * B[k][j];
                }
                C[i][j] = sum;
            }
        }
        return C;
    }

    public static AxisAngle RToAxisAngle(double[][] R) {
        double trace = R[0][0] + R[1][1] + R[2][2];
        double cosTheta = (trace - 1.0) / 2.0;
        cosTheta = Math.max(-1.0, Math.min(1.0, cosTheta));
        double theta = Math.acos(cosTheta);

        double eps = 1e-8;
        double[] axis = new double[3];
        if (theta < eps) {
            axis[0] = 1.0; axis[1] = 0.0; axis[2] = 0.0;
            theta = 0.0;
        } else {
            double denom = 2.0 * Math.sin(theta);
            axis[0] = (R[2][1] - R[1][2]) / denom;
            axis[1] = (R[0][2] - R[2][0]) / denom;
            axis[2] = (R[1][0] - R[0][1]) / denom;
            double n = Math.sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
            axis[0] /= n; axis[1] /= n; axis[2] /= n;
        }
        return new AxisAngle(axis, theta);
    }

    public static class AxisAngle {
        public final double[] axis;
        public final double theta;
        public AxisAngle(double[] axis, double theta) {
            this.axis = axis;
            this.theta = theta;
        }
    }
}
      

For robotics applications using Apache Commons Math, the class org.apache.commons.math3.geometry.euclidean.threed.Rotation provides constructors from axis–angle data and exposes methods to retrieve the axis and angle from a rotation object, implementing the same mathematics as above.

MATLAB provides direct support for axis–angle in the Robotics System Toolbox via axang2rotm and rotm2axang. A from-scratch implementation is:

function R = axisAngleToR(axis, theta)
    % Ensure column vector and normalize
    axis = axis(:);
    axis = axis / norm(axis);

    wx = axis(1); wy = axis(2); wz = axis(3);
    K = [    0, -wz,  wy;
          wz,    0, -wx;
         -wy,  wx,   0];

    R = eye(3) + sin(theta) * K + (1 - cos(theta)) * (K * K);
end

function [axis, theta] = RToAxisAngle(R)
    traceR = trace(R);
    cosTheta = (traceR - 1) / 2;
    cosTheta = max(min(cosTheta, 1), -1);
    theta = acos(cosTheta);
    epsVal = 1e-8;

    if theta < epsVal
        axis = [1; 0; 0];
        theta = 0;
        return;
    end

    axis = 1 / (2 * sin(theta)) * ...
        [R(3,2) - R(2,3);
         R(1,3) - R(3,1);
         R(2,1) - R(1,2)];
    axis = axis / norm(axis);
end
      

In Simulink, one can implement these formulas inside a MATLAB Function block, expose the axis and angle as inputs, and output the rotation matrix R to be used in kinematic chains or visualization blocks. The built-in axang2rotm/rotm2axang functions can also be called directly in such blocks when the Robotics System Toolbox is available.

Mathematica has built-in support for 3D rotations via RotationMatrix and related functions. A symbolic implementation parallel to our derivation is:

skew[{wx_, wy_, wz_}] := {
    {0,   -wz,  wy},
    {wz,   0,  -wx},
    {-wy,  wx,  0}
};

axisAngleToR[axis_List, theta_] := Module[{w, K},
    w = axis/Norm[axis];
    K = skew[w];
    IdentityMatrix[3] + Sin[theta] K + (1 - Cos[theta]) K.K
];

RToAxisAngle[R_] := Module[{trace, cosTheta, theta, axis, eps = 10^-8},
    trace = Tr[R];
    cosTheta = (trace - 1)/2;
    cosTheta = Min[1, Max[-1, cosTheta]];
    theta = ArcCos[cosTheta];

    If[theta < eps,
        Return[{ {1, 0, 0}, 0.0 }]
    ];

    axis = (1/(2 Sin[theta])) * {
        R[[3, 2]] - R[[2, 3]],
        R[[1, 3]] - R[[3, 1]],
        R[[2, 1]] - R[[1, 2]]
    };
    axis = axis/Norm[axis];
    {axis, theta}
];
      

Alternatively, RotationMatrix[theta, axis] constructs the rotation matrix directly, and RotationMatrixToAxisAngle (in newer versions) recovers the axis–angle representation, encapsulating the logarithm on \( \mathrm{SO}(3) \).

This lesson introduced the axis–angle representation as a minimal, geometrically meaningful parameterization of orientations in \( \mathrm{SO}(3) \). We used the skew-symmetric hat operator to connect rotation vectors with the Lie algebra \( \mathfrak{so}(3) \), derived Rodrigues’ formula from the exponential map, and constructed the matrix logarithm to recover axis–angle data from a rotation matrix. We also examined non-uniqueness at \( \theta = \pi \) and small-angle limits, and provided concrete implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica that will serve as building blocks for forward kinematics and more advanced modeling in subsequent chapters.