constructing rotation matrices and generating RELION compatible euler angles from them

aysecan.py

import numpy as np
from scipy.spatial.transform import Rotation as R

# set up rotation matrices for rotations around X axis
rotation_angles = np.linspace(0, 90, 50)
Rx = R.from_euler('x', angles=rotation_angles, degrees=True).as_matrix()
print(Rx.shape)

# set up column vector which points along y
z_vec = np.array([0, 0, 1]).reshape((3, 1))
print(z_vec.shape)

# rotate our z vector by each of our 50 rotations

rotated_z = Rx @ z_vec
print(rotated_z.shape)

# import napari
# viewer = napari.Viewer(ndisplay=3)
# viewer.add_points([0, 0, 0], size=0.1)
# viewer.add_points(rotated_z.reshape((50, 3)))
# napari.run()

# intuition: z vector after rotation is the 3rd column of the rotation matrix... cool!
# if we know what the z vector we want it, we can place it in that third column

# okay, so if we want to construct a matrix, we need three vectors not just one
# the x and y vectors must be orthogonal to both z and each other
# they must also have unit magnitude

# the dot product tells us the projection of a vector onto another vector
# we will initialise a random vector, it might be perpendicular or not
# we will figure out how much of this vector is pointing in the same direction as our plane normal
# we will subtract anything pointing in the direction of the plane normal to get a perpendicular vector

# we have our z vector already
z_vec = np.array([0, 0, 1])

# initialise a random unit vector
random_vector = np.random.random(3)
random_vector /= np.linalg.norm(random_vector)
print(np.linalg.norm(random_vector))

# find the projection of our random vector onto our z vector
projection = np.dot(random_vector, z_vec)
z_component = projection * z_vec
perpendicular_vector = random_vector - z_component
perpendicular_vector /= np.linalg.norm(perpendicular_vector)
print(perpendicular_vector)


# cool! we have two vectors
# almost there

# how do we get vector number 3
# cross product between two vectors gives us a vector perpendicular to both of them
# (perpendicular to the plane containing both of them)

other_perpendicular_vector = np.cross(perpendicular_vector, z_vec)  # order matters here, be careful! right hand rule
other_perpendicular_vector /= np.linalg.norm(other_perpendicular_vector)

# import napari
# viewer = napari.Viewer()
# viewer.add_points(z_vec.reshape(1, 3), size=0.1)
# viewer.add_points(other_perpendicular_vector.reshape((1, 3)), size=0.1, face_color='orange')
# viewer.add_points(perpendicular_vector.reshape((1, 3)), size=0.1, face_color='cornflowerblue')
# napari.run()


# we have constructed a right handed coordinate system, we are great
rotation_matrix = np.zeros((3, 3))
rotation_matrix[:, 0] = other_perpendicular_vector
rotation_matrix[:, 1] = perpendicular_vector
rotation_matrix[:, 2] = z_vec

# we need to calculate RELION format euler angles from this matrix... ugh
# the theory: https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf

# the practice!
euler_angles = R.from_matrix(rotation_matrix).inv().as_euler(seq='ZYZ', degrees=True)

# and because it's useful for you to have
rotation_matrix = R.from_euler(angles=euler_angles, seq='ZYZ', degrees=True).inv().as_matrix()

print(euler_angles)
print(rotation_matrix)