It’s not as simple as you’d think.
If you want to skip the small talk, the code is at the bottom. Sampling 2D rotations uniformly is simple: rotate by an angle from the uniform distribution . Extending this idea to 3D rotations, we could sample each of the three Euler angles from the same uniform distribution . This, however, gives more probability density to transformations which are clustered towards the poles:
In Fast Random Rotation Matrices (James Avro, 1992), a method for uniform random 3D rotation matrices is outlined, the main steps being:
- A random rotation about the z axis
- Rotate the (unmoved) z unit vector to a random position on the unit sphere
The first of these is a simple 2D rotation matrix with a random angle, parameterised by
The second point is more tricky: it uses the Householder matrix with certain properties, given by:
where and are also drawn from . The final, completed random 3D rotation matrix is then (including a reflection through the origin to turn the Householder matrix transformation from a reflection into a rotation):
I have included a python implementation which rotates any vector about the origin with the uniform 3D rotation matrix:
import numpy as np def uniform_random_rotation(x): """Apply a random rotation in 3D, with a distribution uniform over the sphere. Arguments: x: vector or set of vectors with dimension (n, 3), where n is the number of vectors Returns: Array of shape (n, 3) containing the randomly rotated vectors of x, about the mean coordinate of x. Algorithm taken from "Fast Random Rotation Matrices" (James Avro, 1992): https://doi.org/10.1016/B978-0-08-050755-2.50034-8 """ def generate_random_z_axis_rotation(): """Generate random rotation matrix about the z axis.""" R = np.eye(3) x1 = np.random.rand() R[0, 0] = R[1, 1] = np.cos(2 * np.pi * x1) R[0, 1] = -np.sin(2 * np.pi * x1) R[1, 0] = np.sin(2 * np.pi * x1) return R # There are two random variables in [0, 1) here (naming is same as paper) x2 = 2 * np.pi * np.random.rand() x3 = np.random.rand() # Rotation of all points around x axis using matrix R = generate_random_z_axis_rotation() v = np.array([ np.cos(x2) * np.sqrt(x3), np.sin(x2) * np.sqrt(x3), np.sqrt(1 - x3) ]) H = np.eye(3) - (2 * np.outer(v, v)) M = -(H @ R) x = x.reshape((-1, 3)) mean_coord = np.mean(x, axis=0) return ((x - mean_coord) @ M) + mean_coord @ M