English

# Python Script to Calculate Euler Angles and Rotation Matrix

Let’s write a python script to calculate Euler angles and rotation matrix.

## Calculate Rotation Matrix from Euler Angles

This script calculate the 3×3 rotation matrix from Euler angles and rotation order.

=> Table of Euler angles and Rotation matrix

euler.py

import numpy as np

def rotation_matrix(theta1, theta2, theta3, order='xyz'):
"""
input
theta1, theta2, theta3 = rotation angles in rotation order (degrees)
oreder = rotation order of x,y,z　e.g. XZY rotation -- 'xzy'
output
3x3 rotation matrix (numpy array)
"""
c1 = np.cos(theta1 * np.pi / 180)
s1 = np.sin(theta1 * np.pi / 180)
c2 = np.cos(theta2 * np.pi / 180)
s2 = np.sin(theta2 * np.pi / 180)
c3 = np.cos(theta3 * np.pi / 180)
s3 = np.sin(theta3 * np.pi / 180)

if order == 'xzx':
matrix=np.array([[c2, -c3*s2, s2*s3],
[c1*s2, c1*c2*c3-s1*s3, -c3*s1-c1*c2*s3],
[s1*s2, c1*s3+c2*c3*s1, c1*c3-c2*s1*s3]])
elif order=='xyx':
matrix=np.array([[c2, s2*s3, c3*s2],
[s1*s2, c1*c3-c2*s1*s3, -c1*s3-c2*c3*s1],
[-c1*s2, c3*s1+c1*c2*s3, c1*c2*c3-s1*s3]])
elif order=='yxy':
matrix=np.array([[c1*c3-c2*s1*s3, s1*s2, c1*s3+c2*c3*s1],
[s2*s3, c2, -c3*s2],
[-c3*s1-c1*c2*s3, c1*s2, c1*c2*c3-s1*s3]])
elif order=='yzy':
matrix=np.array([[c1*c2*c3-s1*s3, -c1*s2, c3*s1+c1*c2*s3],
[c3*s2, c2, s2*s3],
[-c1*s3-c2*c3*s1, s1*s2, c1*c3-c2*s1*s3]])
elif order=='zyz':
matrix=np.array([[c1*c2*c3-s1*s3, -c3*s1-c1*c2*s3, c1*s2],
[c1*s3+c2*c3*s1, c1*c3-c2*s1*s3, s1*s2],
[-c3*s2, s2*s3, c2]])
elif order=='zxz':
matrix=np.array([[c1*c3-c2*s1*s3, -c1*s3-c2*c3*s1, s1*s2],
[c3*s1+c1*c2*s3, c1*c2*c3-s1*s3, -c1*s2],
[s2*s3, c3*s2, c2]])
elif order=='xyz':
matrix=np.array([[c2*c3, -c2*s3, s2],
[c1*s3+c3*s1*s2, c1*c3-s1*s2*s3, -c2*s1],
[s1*s3-c1*c3*s2, c3*s1+c1*s2*s3, c1*c2]])
elif order=='xzy':
matrix=np.array([[c2*c3, -s2, c2*s3],
[s1*s3+c1*c3*s2, c1*c2, c1*s2*s3-c3*s1],
[c3*s1*s2-c1*s3, c2*s1, c1*c3+s1*s2*s3]])
elif order=='yxz':
matrix=np.array([[c1*c3+s1*s2*s3, c3*s1*s2-c1*s3, c2*s1],
[c2*s3, c2*c3, -s2],
[c1*s2*s3-c3*s1, c1*c3*s2+s1*s3, c1*c2]])
elif order=='yzx':
matrix=np.array([[c1*c2, s1*s3-c1*c3*s2, c3*s1+c1*s2*s3],
[s2, c2*c3, -c2*s3],
[-c2*s1, c1*s3+c3*s1*s2, c1*c3-s1*s2*s3]])
elif order=='zyx':
matrix=np.array([[c1*c2, c1*s2*s3-c3*s1, s1*s3+c1*c3*s2],
[c2*s1, c1*c3+s1*s2*s3, c3*s1*s2-c1*s3],
[-s2, c2*s3, c2*c3]])
elif order=='zxy':
matrix=np.array([[c1*c3-s1*s2*s3, -c2*s1, c1*s3+c3*s1*s2],
[c3*s1+c1*s2*s3, c1*c2, s1*s3-c1*c3*s2],
[-c2*s3, s2, c2*c3]])

return matrix

## Calculate Euler angles from a Rotation Matrix

This script calculate Euler angle from a 3×3 rotation matrix and rotation order.

=> Table of Rotation matrix and Euler angles

euler.py continued

import numpy as np is already done.

def rotation_angles(matrix, order):
"""
input
matrix = 3x3 rotation matrix (numpy array)
oreder(str) = rotation order of x, y, z : e.g, rotation XZY -- 'xzy'
output
theta1, theta2, theta3 = rotation angles in rotation order
"""
r11, r12, r13 = matrix[0]
r21, r22, r23 = matrix[1]
r31, r32, r33 = matrix[2]

if order == 'xzx':
theta1 = np.arctan(r31 / r21)
theta2 = np.arctan(r21 / (r11 * np.cos(theta1)))
theta3 = np.arctan(-r13 / r12)

elif order == 'xyx':
theta1 = np.arctan(-r21 / r31)
theta2 = np.arctan(-r31 / (r11 *np.cos(theta1)))
theta3 = np.arctan(r12 / r13)

elif order == 'yxy':
theta1 = np.arctan(r12 / r32)
theta2 = np.arctan(r32 / (r22 *np.cos(theta1)))
theta3 = np.arctan(-r21 / r23)

elif order == 'yzy':
theta1 = np.arctan(-r32 / r12)
theta2 = np.arctan(-r12 / (r22 *np.cos(theta1)))
theta3 = np.arctan(r23 / r21)

elif order == 'zyz':
theta1 = np.arctan(r23 / r13)
theta2 = np.arctan(r13 / (r33 *np.cos(theta1)))
theta3 = np.arctan(-r32 / r31)

elif order == 'zxz':
theta1 = np.arctan(-r13 / r23)
theta2 = np.arctan(-r23 / (r33 *np.cos(theta1)))
theta3 = np.arctan(r31 / r32)

elif order == 'xzy':
theta1 = np.arctan(r32 / r22)
theta2 = np.arctan(-r12 * np.cos(theta1) / r22)
theta3 = np.arctan(r13 / r11)

elif order == 'xyz':
theta1 = np.arctan(-r23 / r33)
theta2 = np.arctan(r13 * np.cos(theta1) / r33)
theta3 = np.arctan(-r12 / r11)

elif order == 'yxz':
theta1 = np.arctan(r13 / r33)
theta2 = np.arctan(-r23 * np.cos(theta1) / r33)
theta3 = np.arctan(r21 / r22)

elif order == 'yzx':
theta1 = np.arctan(-r31 / r11)
theta2 = np.arctan(r21 * np.cos(theta1) / r11)
theta3 = np.arctan(-r23 / r22)

elif order == 'zyx':
theta1 = np.arctan(r21 / r11)
theta2 = np.arctan(-r31 * np.cos(theta1) / r11)
theta3 = np.arctan(r32 / r33)

elif order == 'zxy':
theta1 = np.arctan(-r12 / r22)
theta2 = np.arctan(r32 * np.cos(theta1) / r22)
theta3 = np.arctan(-r31 / r33)

theta1 = theta1 * 180 / np.pi
theta2 = theta2 * 180 / np.pi
theta3 = theta3 * 180 / np.pi

return (theta1, theta2, theta3)

## Examples

### Calculate Rotation Matrix ⇔ Angles

input in python interpreter



import euler

rotation_mat = euler.rotation_matrix(10, 20, 30, 'yzx')
# calculate a rotation matrix of 10° 20° 30° rotation in YZX order.
print(rotation_mat)
# [[ 0.92541658 -0.20487413  0.31879578]
#  [ 0.34202014  0.81379768 -0.46984631]
#  [-0.16317591  0.54383814  0.82317294]]

angles = euler.rotation_angles(rotation_mat, 'yzx')
#Calculate Euler angles of YZX order from the rotation matrix

print(angles)
# (10.0, 20.0, 29.999999999999993)
#


### Rotate a STL model

A airplane STL is rotated by rotation matrix calculated from Euler angles.
The rotated models are shown in Meshlab.

numpy-stl is used
numpy-stl is a library to manipulate stl and installed by ‘pip’ command.

pip install numpy-stl


import numpy as np
from stl import mesh
# use numpy-stl

airplane = mesh.Mesh.from_file('airplane.stl')

#make a 4x4 matrix and 3x3 rotation matrix is substituted
rotation_x = np.identity(4)
rotation_x[:3, :3] = euler.rotation_matrix(30, 0, 0, 'xzy')

#rotate and save the model
airplane.transform(rotation_x)
airplane.save('airplane_x.stl')
#the model rotated around X axis

#in the same way, rotate the model around Z and Y axis
airplane = mesh.Mesh.from_file('airplane.stl')
# the original stl model needs to be re-read
rotation_xz = np.identity(4)
rotation_xz[:3, :3] = euler.rotation_matrix(30, 60, 0, 'xzy')
airplane.transform(rotation_xz)
airplane.save('airplane_xz.stl')
#the model rotated around X and Z axes

airplane = mesh.Mesh.from_file('airplane.stl')
rotation_xzy = np.identity(4)
rotation_xzy[:3, :3] = euler.rotation_matrix(30, 60, 90, 'xzy')
airplane.transform(rotation_xzy)
airplane.save('airplane_xzy.stl')
#The model rotated around X, Z, and Y axes


local XYZ axes of the airplane

gray： the original model(airplane.stl)
green： the model rotated around X axis(airplane_x.stl)

green： the model rotated around X axis(airplane_x.stl)
orange： the model rotated around X and Z axes(airplane_xz.stl)

orange： the model rotated around X and Z axes(airplane_xz.stl)
yellow： the model rotated around X, Z, and Y axes(airplane_xzy.stl)

## Rotation matrix calculated from three single rotation around XYZ axes

Rotation of the object around the local axes is the same as
Reverse rotation of the global world around the local object

The reverse rotation matrix can be calculated from three single rotations of $$-\theta_1, -\theta_2, -\theta_3$$.

=> Intuition for “rotation around the Local Axis”

euler.py continued

“import numpy as np” is written previously



def axis_rotation(theta, axis_name):
#calculate single rotation of $$\theta$$ matrix around x,y or z
"""
input
theta = rotation angle(degrees)
axis_name = 'x', 'y' or 'z'
output
3x3 rotation matrix
"""

c = np.cos(theta * np.pi / 180)
s = np.sin(theta * np.pi / 180)
if axis_name =='x':
rotation_matrix = np.array([[1, 0, 0],
[0, c, -s],
[0, s, c]])
if axis_name =='y':
rotation_matrix = np.array([[c, 0, s],
[0, 1, 0],
[-s, 0, c]])
elif axis_name =='z':
rotation_matrix = np.array([[c, -s, 0],
[s, c, 0],
[0, 0, 1]])
return rotation_matrix

def confirm_matrix(theta1, theta2, theta3, order='xyz'):
#Calculate Rotation matrix of Euler angles from three single rotation
"""
input
theta1, theta2, theta3 = rotation angles in rotation order
order = rotation order　e.g. XZY rotation -- 'xzy'
output
3x3 rotation matrix
"""
axis_name1 = order[0]
axis_name2 = order[1]
axis_name3 = order[2]

rotation1 = axis_rotation(-theta1, axis_name1)
rotation2 = axis_rotation(-theta2, axis_name2)
rotation3 = axis_rotation(-theta3, axis_name3)

# rotation matrix of global world rotation around the local object
matrix_local = np.dot(rotation3, np.dot(rotation2, rotation1))

# inverse matrix for rotation of the local object
matrix_global = matrix_local.T

return matrix_global

In Python interpreter,

print(euler.rotation_matrix(30, 60, 90, ‘xzy’))
# [[ 3.06161700e-17 -8.66025404e-01 5.00000000e-01]
# [ 5.00000000e-01 4.33012702e-01 7.50000000e-01]
# [-8.66025404e-01 2.50000000e-01 4.33012702e-01]]

print(euler.confirm_matrix(30, 60, 90, ‘xzy’))
# [[ 3.06161700e-17 -8.66025404e-01 5.00000000e-01]
# [ 5.00000000e-01 4.33012702e-01 7.50000000e-01]
# [-8.66025404e-01 2.50000000e-01 4.33012702e-01]]

The output is the same as in the previous script

Rotation matrix can be understood from
animation and numerical formula.
Let’s manipulate the STL model and write numerical formula at the same time for full understanding.