How to change my Code from Euler Angles to Quaternion (or somethin else with uniform distribution)

Question

I have a code, with this i can generate a Volume and put random cubes in there. Now i want these cubes to be randomly rotated. I used euler angles, just because i managed to make it work fast. Now i've read, that euler angles wont produce a really random distributed rotation even if the input is random.

Do you have any idea how i could fix that? Any help is appreciated!

Greetings Paul

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import tifffile as tif


def rotate(Xorg, Yorg, Zorg, alpha, beta,
           gamma):  # TODO:alpha is around y axis, beta is around z axis, gamma is around new z axis
    X_ = Xorg * (np.cos(beta) * np.cos(gamma)) + \
         Yorg * (np.sin(alpha) * np.sin(beta) * np.cos(gamma) - np.cos(alpha) * np.sin(gamma)) + \
         Zorg * (np.cos(alpha) * np.sin(beta) * np.cos(gamma) + np.sin(alpha) * np.sin(gamma))
    # Yaw
    Y_ = Xorg * (np.cos(beta) * np.sin(gamma)) + \
         Yorg * (np.sin(alpha) * np.sin(beta) * np.sin(gamma) + np.cos(alpha) * np.cos(gamma)) + \
         Zorg * (np.cos(alpha) * np.sin(beta) * np.sin(gamma) - np.sin(alpha) * np.cos(gamma))
    # Pitch
    Z_ = Xorg * (-np.sin(beta)) + \
         Yorg * (np.sin(alpha) * np.cos(beta)) + \
         Zorg * (np.cos(alpha) * np.cos(beta))
    # roll
    return X_, Y_, Z_


generate_pic = True
N = 2 ** 8
X = np.linspace(-1, 1, N)
Y = np.linspace(-1, 1, N)
Z = np.linspace(-1, 1, N)
r = 2 ** -4

print("prep done")

X, Y, Z = np.meshgrid(X, Y, Z)

print("meshgrid done")

V = np.zeros((N, N, N), dtype=bool)

for i in range(2**7):
    x0, y0, z0 = 2 * np.random.rand(3) - 1
    alpha, beta, gamma = 2 * np.pi * np.random.rand(3)
    X_, Y_, Z_ = rotate(X, Y, Z, alpha, beta, gamma)
    K = np.logical_and(np.abs(Z_ - z0) < r,
                       np.logical_and(np.abs(X_ - x0) < r, np.abs(Y_ - y0) < r))
    CheckOverlap = np.logical_and(K, V)
    if not CheckOverlap.any():
        V = np.logical_or(K, V)
    else:
        print("Overlap detected")

print("cubes done")
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

z, x, y = V.nonzero()
ax.scatter(x, y, z, c=z, alpha=1)
plt.show()