Intrinsic and Extrinsic Rotations using Euler Angles
I’ve been struggling for a while with 3D rotations, which are very important tool for many applications in physics. In this notebook, we will go on a short tour about Euler angles and in particular, the Bunge convention.
Our goal with this post will be to replicate an animation from Wikipedia:
To me, this an interesting image because it features a rotating coordinate system and a 3D object that changes its orientation based on rotation about several axes.
How can we approach this problem? The ingredients we will need are:
- rotation matrices
- change of basis formulas
- 3D meshing capabilities
Let’s dive in.
Rotation matrices and the change of basis formula¶
Starting with a Z-rotation¶
One of the simplest ways to approach rotations in 3D is through the description of a rotating system around the z-axis.
The starting frame is x-y-z and it becomes, through a rotation of angle $\theta$, the frame x’-y’-z.
A second drawing makes it easy, using trigonometry, to describe the relationship between the rotated basis vectors and the old ones:
Writing the new basis vectors as $e_x’, e_y’, e_z’$, they can be written as functions of the old basis vectors $e_x, e_y, e_z$:
$$ \left \lbrace \begin{aligned} e_x’ = &&\cos \theta e_x + \sin \theta e_y \\ e_y’ = &&-\sin \theta e_x + \cos \theta e_y \\ e_z’ = &&e_z \end{aligned} \right. $$
Change of basis using the R matrix¶
Using this formula, we can ask a first question: given the coordinates of a point $v_x e_x’ + v_y e_y’ + v_z e_z’$ in the rotated frame, what are its coordinates in the original frame?
Plugging these in and applying the above equations, we obtain the coordinates that follow:
$$ v_x e_x’ + v_y e_y’ + v_z e_z’ = u_x e_x + u_y e_y + u_z e_z \\ \leftrightarrow \left \lbrace \begin{aligned} u_x = &&\cos \theta v_x - \sin \theta v_y \\ u_y = &&\sin \theta v_x + \cos \theta v_y \\ u_z = &&v_z \end{aligned} \right. $$
This equation can be reformulated using the familiar rotation matrix $R$ into the equation $u = R v$:
$$ R = \left[ \begin{array}{ccc} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array} \right], \quad u = \left[ \begin{array}{c} u_x \\ u_y \\ u_z \end{array} \right] = R \left[ \begin{array}{c} v_x \\ v_y \\ v_z \end{array} \right] = R v $$
In other words,
the coordinates of the vector $v$ defined in the rotated basis become the coordinates $u = R v$ in the old basis
Inverting the change of basis¶
The relationship can also be inverted, since rotation by amount $\theta$ is the inverse of rotating by $-\theta$. The inverse of a rotation matrix is its transposed matrix, so we can also write:
the coordinates of the vector $u$ defined in the old basis become the coordinates $v = \left[ \begin{array}{ccc} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array} \right]v = R^\top u$ in the rotated basis
It should be noted that although we are talking about rotations, the vectors themselves don’t move, only the frames in which we assign them coordinates. This is usually called a passive rotation, as opposed to the active interpretation of the the same equation. The passive interpretation of the rotation matrix really is an interpretation as a change of coordinate frame, without anything moving.
The active interpretation¶
Where things get hairy for me is that the very same rotation matrix can be used to represent an active rotation, with a moving object but fixed coordinate frames: in that case, the matrix $R$ is a way to compute the changed coordinates of an object after it hase been rotated.
a vector $u’$ defined in the old basis gets rotated and becomes $v’ = R u’$ still in the old basis
This is also interesting, but notice how in that case $u’$ and $v’$ have switched sides when compared with the first formula I gave!
This is definitely a pain point when working with rotation matrices: the matrix doesn’t encode the full information. Knowing on which side the matrix $R$ appears is not enough, you need to know what the vectors on each side stand for, which is either coordinates in different frames (passive interpretation) or coordinates in the same frame (active interpretation).
Putting it all together: Euler angles and the Bunge convention¶
Enough theory! Let’s apply our newfound knowledge. Above, we only discussed a single rotation. In general, we can reach any orientation in 3D space by chaining three rotations around different axes. In my particular field, material science, the way to go is (proper) Euler angles, which are understood to be a triplet of axes around which an object is rotated, two axes being repeated, for example ZXZ.
In particular, the Bunge Euler angle convention used in material science is:
- rotate around axis Z by angle $\varphi_1 \in [0, 2\pi]$
- rotate around axis X’ by angle $\Phi \in [0, \pi]$
- rotate around axis Z’’ by angle $\varphi_2 \in [0, 2\pi]$
Let’ implement this!
In addition to the rotation matrix R around Z, we need a rotation matrix around X, which can be derived easily. Let’s code the matrices using sympy and display them:
import sympy as sp
sp.init_printing()
phi1, Phi, phi2 = sp.symbols('varphi_1 Phi varphi_2')
#phi1, Phi, phi2 = sp.symbols('varphi theta psi')
R_phi1 = sp.Matrix([[sp.cos(phi1), -sp.sin(phi1), 0],
[sp.sin(phi1), sp.cos(phi1), 0],
[0, 0, 1]])
R_Phi = sp.Matrix([[1, 0, 0],
[0, sp.cos(Phi), -sp.sin(Phi)],
[0, sp.sin(Phi), sp.cos(Phi)]])
R_phi2 = sp.Matrix([[sp.cos(phi2), -sp.sin(phi2), 0],
[sp.sin(phi2), sp.cos(phi2), 0],
[0, 0, 1]])
R_phi1, R_Phi, R_phi2
Using the matrices we have just defined, we can formulate a global rotation matrix that, when we multiply a point in 3D space with it, gives us its rotated coordinates in the fixed axis system (notice the primes in the formula below, the same notation as in the theory section):
$$ v’ = R_{\varphi_2} R_{\Phi} R_{\varphi_1} u’ $$
rot_matrix = R_phi2 * R_Phi * R_phi1
rot_matrix
This expression compares favorably with the one in the textbooks. Yay!
Notice also how the order of matrix operations is reversed with respect to the angle triplet $(\varphi_1, \Phi, \varphi_2)$. This is due to the fact that we right-multiply vectors with this matrix to rotate them, see below.
Active rotations in a fixed coordinate frame¶
Let’s now use the matrix above to rotate some objects. As said above, the general formula is:
$$ v’ = R_{\varphi_2} R_{\Phi} R_{\varphi_1} u’, $$
where $u’$ is a vector holding the coordinates of a point and $v’$ are its rotated coordinates, in the same coordinate system. Let’s use pyvista to animate the rotation of an object in a rotation around Z.
import numpy as np
from scipy.spatial.transform import Rotation as R
# convert the rotation matrix so that it can be used with numpy
rot_matrix_numpy = sp.lambdify([phi1, Phi, phi2], rot_matrix)
# double-check our result using scipy
rot_matrix_scipy = lambda phi1, Phi, phi2: R.from_euler('zxz', [phi1, Phi, phi2], degrees=False).as_matrix()
np.allclose(rot_matrix_numpy(1, 2, 3), rot_matrix_scipy(1, 2, 3))
True
from IPython.display import Video
import tempfile
import pyvista as pv
from pyvista import examples
obj = examples.download_teapot().rotate_x(90, inplace=True)
p0 = obj.points.copy()
p = pv.Plotter()
f = tempfile.NamedTemporaryFile(delete=False, suffix='.mp4')
p.open_movie(f.name)
p.add_mesh(obj, style='wireframe')
p.show(interactive=True, auto_close=False, jupyter_backend="none")
for i in range(51):
alpha = 2 * np.pi * i / 50
p.add_title(f"phi_1: {np.rad2deg(alpha):.2f}°")
obj.points = np.array([np.dot(rot_matrix_numpy(alpha, 0, 0), p) for p in p0])
p.write_frame()
p.close()
v = Video(f.name, embed=True)
f.close()
v