raw Math
RAW Math Linear Algebra Affine Transformations

Extracting Translation, Rotation, and Scaling from Homogeneous Transformations

Robert Eisele

A \(4 \times 4\) homogeneous transformation matrix encapsulates translation, rotation, and scaling operations in a single object. It is defined as:

\[\mathbf{T} = \begin{bmatrix} \mathbf{R} & \mathbf{t} \\ \mathbf{0} & 1 \end{bmatrix}\]

where:

Given only \(\mathbf{T}\), how do you recover the translation, the rotation, and the scale that went into it? The translation part is trivial, but the rotation and scaling are entangled inside \(\mathbf{R}\) and need to be pulled apart carefully — and, as it turns out, the obvious way of doing that has a subtle ordering bug that is worth understanding.

Extracting Translation

The translation vector \(\mathbf{t}\) can be read directly from the right column of \(\mathbf{T}\):

\[\mathbf{t} = \begin{bmatrix} \mathbf{T}_{14} \\ \mathbf{T}_{24} \\ \mathbf{T}_{34} \end{bmatrix}\]

Extracting Rotation and Scaling

To separate the rotation and scaling components from the \(3 \times 3\) matrix \(\mathbf{R}\), we first check whether \(\mathbf{R}\) is already a valid rotation matrix, i.e. whether it was built without any scaling at all. Only if it is not do we fall back to a Singular Value Decomposition (SVD).

Validating the Rotation Matrix

A square matrix \(\mathbf{R}\) is a valid rotation matrix if it satisfies:

  1. Orthonormality: the columns (and rows) of \(\mathbf{R}\) are mutually perpendicular and have unit length, which holds iff \(\mathbf{R}^T\) equals \(\mathbf{R}^{-1}\), or \[\mathbf{R}^T \mathbf{R} = \mathbf{I}\]
  2. Determinant: the determinant of \(\mathbf{R}\) is \(+1\), so that it preserves lengths, angles, and orientation (a determinant of \(-1\) would mean \(\mathbf{R}\) also contains a reflection): \[\det(\mathbf{R}) = 1\]

If \(\mathbf{R}\) is not a valid rotation matrix, we decompose it using SVD:

\[\mathbf{R} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T\]

where \(\mathbf{U}\) and \(\mathbf{V}\) are orthogonal and \(\mathbf{\Sigma}\) is diagonal with the (non-negative, descending) singular values of \(\mathbf{R}\) on its diagonal. The closest proper rotation to \(\mathbf{R}\) is then

\[\mathbf{R}_{\text{rot}} = \mathbf{U} \mathbf{V}^T,\]

which is always orthonormal since it is a product of two orthogonal matrices — but it is only a genuine rotation, and not a reflection, as long as \(\det(\mathbf{U}\mathbf{V}^T) = +1\). If the original transform did contain a reflection (a negative scale factor along one axis, say), \(\det(\mathbf{U}\mathbf{V}^T)\) comes out as \(-1\) instead. The standard fix is to flip the sign of the column of \(\mathbf{U}\) belonging to the smallest singular value before forming \(\mathbf{R}_{\text{rot}}\); this yields the closest proper rotation in the Frobenius-norm sense, at the cost of no longer reproducing the reflection exactly — which is expected, since a reflection genuinely cannot be written as a rotation composed with a positive scale.

For the scaling factors, it is tempting to just read them off \(\mathbf{\Sigma}\)'s diagonal directly. This is where the subtlety mentioned earlier comes in: \(\mathbf{\Sigma}\) is always sorted by descending magnitude, and \(\mathbf{V}\) silently carries whatever permutation and sign flips are needed to compensate. If the original per-axis scale factors were not already sorted that way, \(\mathbf{\Sigma}\)'s diagonal no longer lines up with the original \(x\), \(y\), \(z\) axes, and reading it directly reports the right numbers in the wrong order. To recover the scale in the original axes, fold \(\mathbf{\Sigma}\) back through \(\mathbf{V}\) instead:

\[\mathbf{S} = \mathbf{V} \mathbf{\Sigma} \mathbf{V}^T.\]

This is precisely the (symmetric, positive semi-definite) scale factor of the polar decomposition \(\mathbf{R} = \mathbf{R}_{\text{rot}}\mathbf{S}\). Its diagonal gives the correct per-axis scale factors in the original coordinate frame, and its off-diagonal entries — zero for the axis-aligned scaling considered here — would flag any shear if \(\mathbf{R}\) contained some.

Extracting Rotation Angles

To turn \(\mathbf{R}_{\text{rot}}\) into rotation angles around the \(x\), \(y\), and \(z\) axes, we use the Euler angle representation with the ZYX convention, i.e. \(\mathbf{R}_{\text{rot}} = \mathbf{R}_z(\alpha)\mathbf{R}_y(\beta)\mathbf{R}_x(\gamma)\). Given

\[ \mathbf{R}_{\text{rot}} = \begin{bmatrix} r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \end{bmatrix} \]

the angles can be extracted as:

This decomposition is only well-defined away from the gimbal-lock configuration \(\beta = \pm\frac{\pi}{2}\), where the \(x\)- and \(z\)-axis rotations collapse onto the same effective rotation and only their sum or difference can be recovered.

Sidestepping Gimbal Lock with a Unit Quaternion

Gimbal lock is not a flaw in how we read \(\mathbf{R}_{\text{rot}}\) but a property of representing a rotation as three sequential angles at all. If what you need \(\mathbf{R}_{\text{rot}}\) for is interpolating between orientations or storing/transmitting one, rather than reading off three named angles, it is usually better to not extract Euler angles in the first place and convert \(\mathbf{R}_{\text{rot}}\) directly into a unit quaternion \(\hat{\mathbf{q}} = (w, x, y, z)\) instead: a quaternion represents the very same rotation without any coordinate-dependent singularity.

The direct route, \(w=\tfrac12\sqrt{1+\operatorname{tr}(\mathbf{R}_{\text{rot}})}\) followed by

\[(x, y, z) = \frac{1}{4w}\left(r_{32}-r_{23},\ r_{13}-r_{31},\ r_{21}-r_{12}\right),\]

only works away from a \(180°\) rotation, where \(\operatorname{tr}(\mathbf{R}_{\text{rot}})\to -1\) and \(w\to 0\), making that division unstable. The numerically stable version instead branches on whichever of \(\operatorname{tr}(\mathbf{R}_{\text{rot}})\), \(r_{11}\), \(r_{22}\), or \(r_{33}\) is largest, and solves for that component of \(\hat{\mathbf{q}}\) first, before dividing to get the rest — e.g. if \(r_{11}\) is the largest, \(x\) is solved first via \(x=\tfrac12\sqrt{1+r_{11}-r_{22}-r_{33}}\), with \(w\), \(y\), and \(z\) then following by division by \(4x\) instead; the \(r_{22}\)- and \(r_{33}\)-largest cases are the same with the roles of \(x,y,z\) cyclically permuted. Whichever branch is taken, the result should still be renormalized to a unit quaternion, \(\hat{\mathbf{q}} \mathrel{/}= |\hat{\mathbf{q}}|\), to absorb any floating-point drift from \(\mathbf{R}_{\text{rot}}\) not being perfectly orthonormal.

Implementation in Python

Here is a Python example that builds a homogeneous transformation matrix from a known rotation and a (deliberately unsorted) scale, then reconstructs the translation, rotation, scale, and Euler angles from it alone:

import numpy as np

def rotation_matrix_to_euler_angles(R):
    """
    Convert a rotation matrix to Euler angles in the ZYX order.
    :param R: 3x3 rotation matrix
    :return: tuple of Euler angles (yaw, pitch, roll)
    """
    assert R.shape == (3, 3), "Input matrix must be 3x3"
    yaw = np.arctan2(R[1, 0], R[0, 0])
    pitch = np.arcsin(-R[2, 0])
    roll = np.arctan2(R[2, 1], R[2, 2])
    return yaw, pitch, roll

def rotation_matrix_to_quaternion(R):
    """
    Convert a rotation matrix to a unit quaternion (w, x, y, z), using
    whichever of the trace/diagonal entries is largest to avoid dividing
    by something close to zero near a 180-degree rotation.
    """
    tr = R[0, 0] + R[1, 1] + R[2, 2]
    if tr > 0:
        s = np.sqrt(tr + 1.0) * 2
        w = 0.25 * s
        x = (R[2, 1] - R[1, 2]) / s
        y = (R[0, 2] - R[2, 0]) / s
        z = (R[1, 0] - R[0, 1]) / s
    elif R[0, 0] > R[1, 1] and R[0, 0] > R[2, 2]:
        s = np.sqrt(1.0 + R[0, 0] - R[1, 1] - R[2, 2]) * 2
        w = (R[2, 1] - R[1, 2]) / s
        x = 0.25 * s
        y = (R[0, 1] + R[1, 0]) / s
        z = (R[0, 2] + R[2, 0]) / s
    elif R[1, 1] > R[2, 2]:
        s = np.sqrt(1.0 + R[1, 1] - R[0, 0] - R[2, 2]) * 2
        w = (R[0, 2] - R[2, 0]) / s
        x = (R[0, 1] + R[1, 0]) / s
        y = 0.25 * s
        z = (R[1, 2] + R[2, 1]) / s
    else:
        s = np.sqrt(1.0 + R[2, 2] - R[0, 0] - R[1, 1]) * 2
        w = (R[1, 0] - R[0, 1]) / s
        x = (R[0, 2] + R[2, 0]) / s
        y = (R[1, 2] + R[2, 1]) / s
        z = 0.25 * s
    q = np.array([w, x, y, z])
    return q / np.linalg.norm(q)  # renormalize to a unit quaternion

# Define the rotation angle in radians
theta = np.radians(30)

# Create a rotation matrix around the Z-axis
R_z = np.array([
    [np.cos(theta), -np.sin(theta), 0],
    [np.sin(theta), np.cos(theta), 0],
    [0, 0, 1]
])

# Introduce scaling factors for each axis, deliberately NOT sorted by magnitude
scaling_factors = np.array([1.2, 0.8, 1.0])
S_diag = np.diag(scaling_factors)

# Combine rotation and scaling
R_combined = np.dot(R_z, S_diag)

# Define a homogeneous transformation matrix with translation
T = np.array([
    [R_combined[0, 0], R_combined[0, 1], R_combined[0, 2], 1],
    [R_combined[1, 0], R_combined[1, 1], R_combined[1, 2], 2],
    [R_combined[2, 0], R_combined[2, 1], R_combined[2, 2], 3],
    [0, 0, 0, 1]
])

## Now try to reconstruct the inputs:

# Extract the rotation and scaling matrix R from T
R = T[:3, :3]

# Check if R already is a valid rotation matrix
if np.allclose(np.dot(R.T, R), np.eye(3)) and np.isclose(np.linalg.det(R), 1.0):
    R_rot = R
    S = np.eye(3)
else:
    # Perform SVD on the extracted matrix
    U, Sigma, Vt = np.linalg.svd(R)

    # Guard against a reflection: force det(U @ Vt) = +1 by flipping the
    # column of U that belongs to the smallest singular value
    if np.linalg.det(U @ Vt) < 0:
        U[:, -1] *= -1

    R_rot = np.dot(U, Vt)

    # Fold Sigma back through V to realign it with the ORIGINAL axes --
    # Sigma alone is sorted by magnitude and not axis-aligned in general
    S = Vt.T @ np.diag(Sigma) @ Vt

# Extract the translation vector t from T
t = T[:3, 3]

# Extract Euler angles from the rotation matrix
yaw, pitch, roll = rotation_matrix_to_euler_angles(R_rot)

# Extract a unit quaternion from the rotation matrix as a gimbal-lock-free
# alternative to the Euler angles above
quat = rotation_matrix_to_quaternion(R_rot)

# Print results
print("Extracted Rotation Matrix R_rot:")
print(R_rot)
print("Extracted Scaling Factors:")
print(np.diag(S))
print("Extracted Translation Vector t:")
print(t)
print("Extracted Euler Angles (yaw, pitch, roll):")
print(np.degrees([yaw, pitch, roll]))
print("Extracted Unit Quaternion (w, x, y, z):")
print(quat)

Running this prints back the exact input scale \([1.2,\ 0.8,\ 1.0]\) in the correct axis order and a \(30°\) yaw with zero pitch and roll — reading \(\mathbf{\Sigma}\)'s diagonal directly instead would instead have reported \([1.2,\ 1.0,\ 0.8]\), silently swapping the \(y\)- and \(z\)-scale. The quaternion comes out as \((\cos 15°,\ 0,\ 0,\ \sin 15°)\), the expected half-angle representation of the same \(30°\) rotation around the \(z\)-axis.

Related Stuff

raw / Geometry
24 Aug
·
Robert Eisele

Introduction to Rotations with Complex Numbers and Quaternions

A rotation is defined by what it preserves, and this chapter follows that definition from the plane into three dimensions. It derives the 2D rotation matrix twice, rewrites it as multiplication by a unit complex number, and shows that composing two reflections produces a rotation by twice the angle between the mirrors, which is where the half angle comes from. Euler's rotation theorem and Frobenius's theorem then explain why three dimensions need a four-dimensional, non-commutative algebra. Rodrigues' rotation formula is derived geometrically and reproduced exactly by the quaternion sandwich product, whose two-sided form, half angle and two-to-one sign ambiguity all follow from the reflection picture. The chapter closes with composition, matrix conversion in both directions including Shepperd's stable branch selection, the efficient two-cross-product evaluation, spherical linear interpolation, angular velocity and integration, and why gimbal lock is a property of Euler angles rather than of rotation. Nine visualizations accompany the derivations, together with JavaScript code using Quaternion.js.
raw / Linear Algebra
02 Aug
·
Robert Eisele

Introduction to Singular Value Decomposition

This derives the singular value decomposition (SVD) of a real matrix of any shape directly from the eigendecomposition of the symmetric positive semidefinite matrix A^T A, motivated visually by how any matrix maps the unit circle to an ellipse and by the general rotate-stretch-rotate geometric picture this gives for every matrix. It works through a full numerical 2x2 example, shows how the right and left singular vectors are eigenvectors of A^T A and A A^T and how this completes the four fundamental subspaces picture, illustrated with a diagram of the Strang diagram itself, distinguishes the full and compact SVD, answers precisely when the eigendecomposition of a symmetric matrix coincides with its SVD (positive semidefiniteness) and how singular values differ from eigenvalues in general (including for defective and rotation matrices), derives the uniqueness of the SVD up to sign flips and rotations within repeated singular values with a worked cautionary example of what goes wrong when both sides are eigendecomposed independently, explains the numerically robust way to derive one side from the other, proves the Eckart-Young theorem for optimal low-rank approximation, and lists its properties (including the determinant as a product of singular values) and applications in least squares (including the explicit pseudo-inverse sum formula, exact-solvability criterion, and affine solution space), homogeneous least squares, data compression, and principal component analysis.
raw / Algebra
21 Feb
·
Robert Eisele

Introduction to Quaternions

Hamiltonian quaternions extend the complex numbers into a non-commutative algebra built from the imaginary units i, j, and k. The treatment develops quaternion notation and arithmetic, conjugation, norm, inverse, division, exponential, principal logarithm, polar decomposition, De Moivre's formula, and scalar and quaternion powers with their ordering and branch conventions. Unit quaternions are then used to represent, compose, and invert 3D rotations without gimbal lock.
All Math