Composition of finite rotations

01

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 $\htmlClass{fea-math-bold-italic}{\boldsymbol l}$: $\htmlClass{fea-math-bold-italic}{\boldsymbol a}=(\htmlClass{fea-math-bold-italic}{\boldsymbol l},\alpha)$, and
• rotation (b) of angle $\beta$ around the unit vector $\htmlClass{fea-math-bold-italic}{\boldsymbol m}$: $\htmlClass{fea-math-bold-italic}{\boldsymbol b}=(\htmlClass{fea-math-bold-italic}{\boldsymbol m},\beta)$.

The composition of these two rotations gives the rotation of angle $\gamma$ around the unit vector $\htmlClass{fea-math-bold-italic}{\boldsymbol n}$: $\htmlClass{fea-math-bold-italic}{\boldsymbol c}=(\htmlClass{fea-math-bold-italic}{\boldsymbol n},\gamma)$.

The result of composition:

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 $\htmlClass{fea-math-bold-italic}{\boldsymbol n}$ are:

where $\htmlClass{fea-math-bold-italic}{\boldsymbol l}\cdot\htmlClass{fea-math-bold-italic}{\boldsymbol m}$ denotes the dot product and $\htmlClass{fea-math-bold-italic}{\boldsymbol l}\times\htmlClass{fea-math-bold-italic}{\boldsymbol m}$ denotes the cross product.¹⁾ Obviously:

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

where the vectors $\htmlClass{fea-math-bold-italic}{\boldsymbol a}_l$ and $\htmlClass{fea-math-bold-italic}{\boldsymbol b}_m$ are:

and:

02

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 $\htmlClass{fea-math-bold-italic}{\boldsymbol a}$ and $\htmlClass{fea-math-bold-italic}{\boldsymbol b}$ as:

and:

Evidently, when a quaternion represents a rotation:

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

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

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

Therefore:

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

The conjugate of:

is:

and:

03

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 $\htmlClass{fea-math-bold-italic}{\boldsymbol R}$ (see Section 21.1). The composition of two rotations $\htmlClass{fea-math-bold-italic}{\boldsymbol R}_1$ and $\htmlClass{fea-math-bold-italic}{\boldsymbol R}_2$ ($\htmlClass{fea-math-bold-italic}{\boldsymbol R}_1$ is applied first) is:

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

Here $\htmlClass{fea-math-bold-italic}{\boldsymbol{\rho}}$ is a three-component rotation vector. The axis-angle parameter form is:

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:

Quaternion representation is one of the most utilized:

where:

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

The equivalent rotation matrix expressed in quaternion elements is:

Obviously, substituting the expressions of $q_0$, $q_1$, $q_2$ and $q_3$ in $\htmlClass{fea-math-bold-italic}{\boldsymbol R}_q$, the $\htmlClass{fea-math-bold-italic}{\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: $\htmlClass{fea-math-bold-italic}{\boldsymbol R}_{\mathit{TB}}=\htmlClass{fea-math-bold-italic}{\boldsymbol R}_{\mathit{active}}^T$. The rotation matrix is:

04

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.

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.

05

Example

Let's consider two rotations:

a) $\alpha=\dfrac{\pi}{4}$ around the direction defined by the vector $\htmlClass{fea-math-bold-italic}{\boldsymbol d}_1=\begin{bmatrix}1&3&2\end{bmatrix}$. Direction cosines: $\htmlClass{fea-math-bold-italic}{\boldsymbol l}=\htmlClass{fea-math-bold-italic}{\boldsymbol d}_1/\|\htmlClass{fea-math-bold-italic}{\boldsymbol d}_1\|$.

b) $\beta=\dfrac{\pi}{3}$ around the direction defined by the vector $\htmlClass{fea-math-bold-italic}{\boldsymbol d}_2=\begin{bmatrix}2&-1&1\end{bmatrix}$. Direction cosines: $\htmlClass{fea-math-bold-italic}{\boldsymbol m}=\htmlClass{fea-math-bold-italic}{\boldsymbol d}_2/\|\htmlClass{fea-math-bold-italic}{\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
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)

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:

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:

Using quaternions (run comp2):

% 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)

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$:

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.

06

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.