21.3. Composition of finite rotations

Composition of finite rotations

Unless otherwise stated, the active rotation convention introduced in Section 21.1 is used throughout this section.

Euler's rotation theorem states: “in three-dimensional space, any displacement of a rigid body such that a point on the rigid body remains fixed, is equivalent to a single rotation about some axis that runs through the fixed point”, [1]. It also follows that the composition of two rotations is a rotation.

In the middle of the nineteenth century (1840) the French mathematician Olinde Rodrigues demonstrated the formula for composition of large rotations, [2]. Let’s consider two rotations:

rotation (a) of angle $\alpha$ around the unit vector $\boldsymbol l$ : $\boldsymbol a=(\boldsymbol l,\alpha)$ , and
rotation (b) of angle $\beta$ around the unit vector $\boldsymbol m$ : $\boldsymbol b=(\boldsymbol m,\beta)$ .

The composition of these two rotations gives the rotation of angle $\gamma$ around the unit vector $\boldsymbol n$ : $\boldsymbol c=(\boldsymbol n,\gamma)$ .

The result of composition:

\boldsymbol c=\boldsymbol a\circ\boldsymbol b\qquad\text{or}\qquad(\boldsymbol n,\gamma)=(\boldsymbol l,\alpha)\circ(\boldsymbol m,\beta)

The composition of finite rotations is not commutative; the result depends on the order of application of the rotations! In the composition considered here, rotation (b) is applied first and rotation (a) second. The angle $\gamma$ and the unit vector $\boldsymbol n$ are:

\begin{aligned}\sin\frac{\gamma}{2}\boldsymbol n&=\sin\frac{\alpha}{2}\cos\frac{\beta}{2}\boldsymbol l+\cos\frac{\alpha}{2}\sin\frac{\beta}{2}\boldsymbol m+\sin\frac{\alpha}{2}\sin\frac{\beta}{2}\boldsymbol l\times\boldsymbol m\\ \cos\frac{\gamma}{2}&=\cos\frac{\alpha}{2}\cos\frac{\beta}{2}-\sin\frac{\alpha}{2}\sin\frac{\beta}{2}\boldsymbol l\cdot\boldsymbol m \end{aligned}\tag{*}

where $\boldsymbol l\cdot\boldsymbol m$ denotes the dot product and $\boldsymbol l\times\boldsymbol m$ denotes the cross product.¹⁾ Obviously:

\boldsymbol l\cdot\boldsymbol l=1\qquad\boldsymbol m\cdot\boldsymbol m=1\qquad\boldsymbol n\cdot\boldsymbol n=1

The two equations (*) may be written in a nicer form:

\boldsymbol c_n=\left(b_0\boldsymbol a_l+a_0\boldsymbol b_m+\boldsymbol a_l\times\boldsymbol b_m\quad a_0b_0-\boldsymbol a_l\cdot\boldsymbol b_m\right)

where the vectors $\boldsymbol a_l$ and $\boldsymbol b_m$ are:

\boldsymbol a_l=\boldsymbol l\sin\frac{\alpha}{2}=\begin{Bmatrix}l_x\\l_y\\l_z\end{Bmatrix}\sin\frac{\alpha}{2}\qquad \boldsymbol b_m=\boldsymbol m\sin\frac{\beta}{2}=\begin{Bmatrix}m_x\\m_y\\m_z\end{Bmatrix}\sin\frac{\beta}{2}\qquad \boldsymbol c_n=\boldsymbol n\sin\frac{\gamma}{2}=\begin{Bmatrix}n_x\\n_y\\n_z\end{Bmatrix}\sin\frac{\gamma}{2}

and:

a_0=\cos\frac{\alpha}{2}\qquad b_0=\cos\frac{\beta}{2}

Quaternions

More elegant expressions are obtained if we represent rotations as quaternions. More about quaternions in [3], [4]. Throughout this section, quaternions are written in scalar-first form. Let’s define two quaternions $\boldsymbol a$ and $\boldsymbol b$ as:

\boldsymbol a=\begin{Bmatrix}a_0\\\boldsymbol a_l\end{Bmatrix}=\begin{Bmatrix}\cos\frac{\alpha}{2}\\\boldsymbol l\sin\frac{\alpha}{2}\end{Bmatrix}=\begin{Bmatrix}\cos\frac{\alpha}{2}\\l_x\sin\frac{\alpha}{2}\\l_y\sin\frac{\alpha}{2}\\l_z\sin\frac{\alpha}{2}\end{Bmatrix}=\begin{Bmatrix}a_0\\a_1\\a_2\\a_3\end{Bmatrix}

and:

\boldsymbol b=\begin{Bmatrix}b_0\\\boldsymbol b_m\end{Bmatrix}=\begin{Bmatrix}\cos\frac{\beta}{2}\\\boldsymbol m\sin\frac{\beta}{2}\end{Bmatrix}=\begin{Bmatrix}\cos\frac{\beta}{2}\\m_x\sin\frac{\beta}{2}\\m_y\sin\frac{\beta}{2}\\m_z\sin\frac{\beta}{2}\end{Bmatrix}=\begin{Bmatrix}b_0\\b_1\\b_2\\b_3\end{Bmatrix}

Evidently, when a quaternion represents a rotation:

\|\boldsymbol a\|=\boldsymbol a^T\boldsymbol a=1\qquad\|\boldsymbol b\|=\boldsymbol b^T\boldsymbol b=1

Quaternion multiplication used here is the Hamilton product. The composition (multiplication) of two quaternions gives, [3], [4]:

\boldsymbol c=\begin{Bmatrix}c_0\\c_1\\c_2\\c_3\end{Bmatrix}=\boldsymbol a\circ\boldsymbol b =\begin{Bmatrix}\cos\frac{\gamma}{2}\\n_x\sin\frac{\gamma}{2}\\n_y\sin\frac{\gamma}{2}\\n_z\sin\frac{\gamma}{2}\end{Bmatrix} =\begin{Bmatrix} a_0b_0-a_1b_1-a_2b_2-a_3b_3\\ a_0b_1+a_1b_0+a_2b_3-a_3b_2\\ a_0b_2-a_1b_3+a_2b_0+a_3b_1\\ a_0b_3+a_1b_2-a_2b_1+a_3b_0 \end{Bmatrix}\tag{**}

Quaternions may also be considered as an extension of complex numbers and they can be written like this:

a=a_0+a_1\boldsymbol i+a_2\boldsymbol j+a_3\boldsymbol k

Composition of two rotations is equivalent to the multiplication of two quaternions written in this form and considering these basic multiplication rules:

\begin{gathered} \boldsymbol i^2=\boldsymbol j^2=\boldsymbol k^2=-1\\ \boldsymbol i\boldsymbol j=\boldsymbol k\qquad\boldsymbol j\boldsymbol i=-\boldsymbol k\\ \boldsymbol j\boldsymbol k=\boldsymbol i\qquad\boldsymbol k\boldsymbol j=-\boldsymbol i\\ \boldsymbol k\boldsymbol i=\boldsymbol j\qquad\boldsymbol i\boldsymbol k=-\boldsymbol j \end{gathered}

Therefore:

\begin{aligned} a\circ b&=(a_0+a_1\boldsymbol i+a_2\boldsymbol j+a_3\boldsymbol k)(b_0+b_1\boldsymbol i+b_2\boldsymbol j+b_3\boldsymbol k)\\ &=a_0b_0+a_0b_1\boldsymbol i+a_0b_2\boldsymbol j+a_0b_3\boldsymbol k\\ &\quad+a_1b_0\boldsymbol i-a_1b_1+a_1b_2\boldsymbol k-a_1b_3\boldsymbol j\\ &\quad+a_2b_0\boldsymbol j-a_2b_1\boldsymbol k-a_2b_2+a_2b_3\boldsymbol i\\ &\quad+a_3b_0\boldsymbol k+a_3b_1\boldsymbol j-a_3b_2\boldsymbol i-a_3b_3\\ &=a_0b_0-a_1b_1-a_2b_2-a_3b_3\\ &\quad+(a_0b_1+a_1b_0+a_2b_3-a_3b_2)\boldsymbol i\\ &\quad+(a_0b_2-a_1b_3+a_2b_0+a_3b_1)\boldsymbol j\\ &\quad+(a_0b_3+a_1b_2-a_2b_1+a_3b_0)\boldsymbol k \end{aligned}

This result is identical to (**). The product of two quaternions is not commutative!

The conjugate of:

a=a_0+a_1\boldsymbol i+a_2\boldsymbol j+a_3\boldsymbol k

is:

a^*=a_0-a_1\boldsymbol i-a_2\boldsymbol j-a_3\boldsymbol k

and:

aa^*=|a|^2=a_0^2+a_1^2+a_2^2+a_3^2

Rotation representations

There are several ways to describe a rotation; four of them are considered here.

The most common and most used is the rotation matrix $\boldsymbol R$ (see Section 21.1). The composition of two rotations $\boldsymbol R_1$ and $\boldsymbol R_2$ ( $\boldsymbol R_1$ is applied first) is:

\boldsymbol R=\boldsymbol R_2\boldsymbol R_1

The axis-angle representation describes a rotation of angle $\varphi$ and direction defined by the unit vector $\boldsymbol n=(n_x\ n_y\ n_z)^T$ . Using the notation introduced in Section 21.1, the vector form is:

\boldsymbol\rho=\varphi\boldsymbol n,\qquad\|\boldsymbol\rho\|=\varphi.

Here $\boldsymbol\rho$ is a three-component rotation vector. The axis-angle parameter form is:

\boldsymbol r=\begin{bmatrix}n_x&n_y&n_z&\varphi\end{bmatrix}.

This is a four-parameter axis-angle array, as used by MATLAB functions such as axang2rotm and rotm2axang. The two forms contain the same physical information but are not the same mathematical object.

The equivalent rotation matrix (see Section 21.1) is expressed in axis-angle parameters:

\boldsymbol R_{aa}=\begin{bmatrix} \cos\varphi+n_x^2(1-\cos\varphi)&n_xn_y(1-\cos\varphi)-n_z\sin\varphi&n_zn_x(1-\cos\varphi)+n_y\sin\varphi\\ n_xn_y(1-\cos\varphi)+n_z\sin\varphi&\cos\varphi+n_y^2(1-\cos\varphi)&n_yn_z(1-\cos\varphi)-n_x\sin\varphi\\ n_zn_x(1-\cos\varphi)-n_y\sin\varphi&n_yn_z(1-\cos\varphi)+n_x\sin\varphi&\cos\varphi+n_z^2(1-\cos\varphi) \end{bmatrix}

Quaternion representation is one of the most utilized:

\boldsymbol q=\begin{Bmatrix}q_0\\q_1\\q_2\\q_3\end{Bmatrix} =\begin{Bmatrix}\cos\frac{\varphi}{2}\\n_x\sin\frac{\varphi}{2}\\n_y\sin\frac{\varphi}{2}\\n_z\sin\frac{\varphi}{2}\end{Bmatrix} =\begin{Bmatrix}\cos\frac{\varphi}{2}\\\boldsymbol n\sin\frac{\varphi}{2}\end{Bmatrix}

where:

\boldsymbol q^T\boldsymbol q=q_0^2+q_1^2+q_2^2+q_3^2=1

A unit quaternion $q$ and its opposite $-q$ represent the same finite rotation.

The equivalent rotation matrix expressed in quaternion elements is:

\boldsymbol R_q=\begin{bmatrix} 1-2(q_2^2+q_3^2)&2(q_1q_2-q_0q_3)&2(q_1q_3+q_0q_2)\\ 2(q_1q_2+q_0q_3)&1-2(q_3^2+q_1^2)&2(q_2q_3-q_0q_1)\\ 2(q_1q_3-q_0q_2)&2(q_2q_3+q_0q_1)&1-2(q_1^2+q_2^2) \end{bmatrix}

Obviously, substituting the expressions of $q_0$ , $q_1$ , $q_2$ and $q_3$ in $\boldsymbol R_q$ , the $\boldsymbol R_{aa}$ rotation matrix is obtained immediately.

Tait–Bryan angles. The three angles, called yaw, pitch and roll, are $(\psi\ \theta\ \varphi)$ , see Section 21.2: $\boldsymbol R_{\mathrm{TB}}=\boldsymbol R_{\mathrm{active}}^T$ . The rotation matrix is:

\boldsymbol R_{\mathrm{TB}}=\begin{bmatrix} \cos\psi\cos\theta&\sin\psi\cos\theta&-\sin\theta\\ -\sin\psi\cos\varphi+\cos\psi\sin\theta\sin\varphi&\cos\psi\cos\varphi+\sin\psi\sin\theta\sin\varphi&\cos\theta\sin\varphi\\ \sin\psi\sin\varphi+\cos\psi\sin\theta\cos\varphi&-\cos\psi\sin\varphi+\sin\psi\sin\theta\cos\varphi&\cos\theta\cos\varphi \end{bmatrix}

MATLAB conversions

MATLAB commands to convert different rotation representations are listed below. For the relationship between quaternions and Tait–Bryan angles, see also [5].

The MATLAB functions listed below use the active rotation convention adopted in Section 21.1. For Euler-angle conversions, MATLAB uses the intrinsic “ZYX” sequence by default, corresponding to yaw–pitch–roll. This is the same axis order used in Section 21.2; however, because Section 21.2 uses the global-to-local matrix convention, the corresponding rotation matrix is the transpose of MATLAB’s active rotation matrix.

From / To	Rotation matrix $R$	Axis-angle $r$	Quaternion $q$	Tait–Bryan angles $e$
Rotation matrix $R$	—	`rotm2axang(R)`	`rotm2quat(R)`	`rotm2eul(R,"ZYX")`
Axis-angle $r$	`axang2rotm(r)`	—	`axang2quat(r)`	—
Quaternion $q$	`quat2rotm(q)`	`quat2axang(q)`	—	`quat2eul(q,"ZYX")`
Tait–Bryan angles $e$	`eul2rotm(e,"ZYX")`	—	`eul2quat(e,"ZYX")`	—

\boldsymbol R_{\mathrm{TB}}=\operatorname{eul2rotm}(e,\text{"ZYX"})^T

e=\operatorname{rotm2eul}(\boldsymbol R_{\mathrm{TB}}^T,\text{"ZYX"})

These modern conversion functions are provided by Robotics System Toolbox.

The legacy functions vrrotmat2vec and vrrotvec2mat are deprecated and scheduled for removal. The modern functions rotm2axang and axang2rotm should be used instead.

Example

Let's consider two rotations:

a) $\alpha=\dfrac{\pi}{4}$ around the direction defined by the vector $\boldsymbol d_1=\begin{bmatrix}1&3&2\end{bmatrix}$ . Direction cosines: $\boldsymbol l=\boldsymbol d_1/\|\boldsymbol d_1\|$ .

b) $\beta=\dfrac{\pi}{3}$ around the direction defined by the vector $\boldsymbol d_2=\begin{bmatrix}2&-1&1\end{bmatrix}$ . Direction cosines: $\boldsymbol m=\boldsymbol d_2/\|\boldsymbol d_2\|$ .

The composition of these two rotations, first rotation (a) and then rotation (b), can be performed with these MATLAB programs:

Using rotation matrices (run comp1):

comp1.m

% comp1
al=pi/4;
l=[1 3 2]; l=l/norm(l);
ra=[l al];

be=pi/3;
m=[2 -1 1]; m=m/norm(m);
rb=[m be];

Ra=axa2R(ra);      % convert axis-angle to rotation matrix
Rb=axa2R(rb);
%Ra=axang2rotm(ra);      % convert axis-angle to rotation matrix (Matlab command)
%Rb=axang2rotm(rb);

R=Rb*Ra;
rc=R2axa(R);       % convert rotation matrix to axis-angle 
%rc=rotm2axang(R); % convert rotation matrix to axis-angle (Matlab command)
n=rc(1:3)'
ga=rc(4)

axa2R.m

function R=axa2R(r)
n13=norm(r(1:3));
r(1:3)=r(1:3)/n13;
nx=r(1); ny=r(2); nz=r(3);
fi=r(4);
R=[  cos(fi)-nx^2*(cos(fi)-1),       - nz*sin(fi) - nx*ny*(cos(fi) - 1), ny*sin(fi) - nx*nz*(cos(fi) - 1)
     nz*sin(fi) - nx*ny*(cos(fi) - 1), cos(fi)-ny^2*(cos(fi)-1),       - nx*sin(fi) - ny*nz*(cos(fi) - 1)
   - ny*sin(fi) - nx*nz*(cos(fi) - 1), nx*sin(fi) - ny*nz*(cos(fi) - 1), cos(fi)-nz^2*(cos(fi)-1)        ];

The function R2axa treats $\gamma\approx0$ , $0<\gamma<\pi$ and $\gamma\approx\pi$ separately, clips the argument of acos to $[-1,1]$ , and normalizes the recovered axis:

R2axa.m

function r=R2axa(R)
%*** R2axa ***
% Convert an active rotation matrix to principal axis-angle form.
tol=1e-12;
c=(trace(R)-1)/2;
c=max(-1,min(1,c));
ga=acos(c);
w=[R(3,2)-R(2,3) R(1,3)-R(3,1) R(2,1)-R(1,2)]/2;
sn=norm(w);
if ga<1e-7
    if sn<tol
        r=[1 0 0 0];
        return
    end
    ga=atan2(sn,c);
    n=w/sn;
elseif pi-ga<1e-5
    % Recover the axis without division by sin(ga).
    A=(R+R'+2*eye(3))/4;
    [~,k]=max(diag(A));
    n=zeros(1,3);
    n(k)=sqrt(max(A(k,k),0));
    j=setdiff(1:3,k);
    n(j)=A(j,k)'/n(k);
    if dot(n,w)<0; n=-n; end
    ga=atan2(sn,c);
else
    n=w/sin(ga);
end
n=n/norm(n);
r=[n ga];
end

The result is: $rc=\begin{bmatrix}0.5766&0.0231&0.8167&1.3547\end{bmatrix}$ , that is:

\boldsymbol n=\begin{Bmatrix}0.5766\\0.0231\\0.8167\end{Bmatrix}\qquad\text{and}\qquad\gamma=1.3547\ \mathrm{rad}

Using quaternions (run comp2):

comp2.m

% comp2
al=pi/4;
l=[1 3 2]; l=l/norm(l);
qa=[cos(al/2) sin(al/2)*l];     % first quaternion 
 
be=pi/3;
m=[2 -1 1]; m=m/norm(m);
qb=[cos(be/2) sin(be/2)*m];     % second quaternion

qc=quatmult(qb,qa);             % quaternion multiplication
%qc=quatmultiply(qb,qa);        % quaternion multiplication (Matlab command)
rc=quat2axa(qc);                % convert quaternion to axis-angle 
%rc=quat2axang(qc);             % convert quaternion to axis-angle (Matlab command)

n=rc(1:3)' 
ga=rc(4)

quatmult.m

function qc=quatmult(qa,qb)
a0=qa(1); a1=qa(2); a2=qa(3); a3=qa(4);
b0=qb(1); b1=qb(2); b2=qb(3); b3=qb(4);
qc=[a0*b0-a1*b1-a2*b2-a3*b3
    a0*b1+a1*b0+a2*b3-a3*b2
    a0*b2-a1*b3+a2*b0+a3*b1
    a0*b3+a1*b2-a2*b1+a3*b0];

The function quat2axa normalizes the quaternion and selects the sign with $q_0\ge0$ . It protects the argument of acos, treats zero rotation explicitly, and returns the principal angle $0\le\gamma\le\pi$ :

quat2axa.m

function r=quat2axa(q)
%*** quat2axa ***
% Convert a scalar-first quaternion to principal axis-angle form.
tol=1e-12;
q=q(:)';
nq=norm(q);
if nq==0; error('The quaternion must be nonzero.'); end
q=q/nq;
if q(1)<0; q=-q; end
c=max(-1,min(1,q(1)));
ga=2*acos(c);
sn=norm(q(2:4));
if sn<tol
    r=[1 0 0 0];
    return
end
if ga<1e-7; ga=2*atan2(sn,c); end
n=q(2:4)/sn;
n=n/norm(n);
r=[n ga];
end

For $\gamma=0$ , the rotation axis is indeterminate; any unit vector represents the identity rotation. In the program, the conventional choice $[1\ 0\ 0]$ is used.

References

[1] Euler's rotation theorem.

[2] Composition of two axis-angle rotations.

[3] Quaternion.

[4] Quaternions and spatial rotation.

[5] Getting Euler–Tait–Bryan angles from quaternion representation.

¹⁾ Dot product and cross product of two vectors in matrix form:

\boldsymbol a=\begin{Bmatrix}a_x\\a_y\\a_z\end{Bmatrix}\qquad\text{and}\qquad\boldsymbol b=\begin{Bmatrix}b_x\\b_y\\b_z\end{Bmatrix}

\boldsymbol a\cdot\boldsymbol b=a_xb_x+a_yb_y+a_zb_z\qquad \boldsymbol a\times\boldsymbol b=\begin{bmatrix}0&-a_z&a_y\\a_z&0&-a_x\\-a_y&a_x&0\end{bmatrix}\begin{Bmatrix}b_x\\b_y\\b_z\end{Bmatrix}