My goal is to uniformly sample the distribution of all unit quaternions $q$ that have angular distance $\theta$ relative to a given quaternion orientation $q_1 = [s_1, q_{1,1}, q_{1,2},q_{1,3}]$.
I calculate angular distance between two unit quaternions $q_1$ and $q_2$ as follows (this is from this post):
$$ \theta = \arccos(2\langle q_1,q_2 \rangle^2 - 1) \tag{1}$$ where $$ \langle q_1,q_2 \rangle = s_1s_2+q_{1,1}q_{2,1}+q_{1,2}q_{2,2}+q_{1,3}q_{2,3} \tag{2} $$
In addition, I know from here that all random unit quaternions describing 3D orientations can be addressed by uniformly sampling $u_1,u_2,u_3 \in [0,1]$ and computing the quaternion $$ q_x = (\sqrt{1-u_1}\sin2\pi u_2,\sqrt{1-u_1}\cos2\pi u_2,\sqrt{u_1}\sin2\pi u_3,\sqrt{u_1}\cos2\pi u_3) \tag{3}$$
I can combine these equations algebraically to get
$$ \sqrt{\frac{\cos(\theta)+1}{2}} = \langle q_1,q_x \rangle \tag{4}$$
Where $q_x$ is as defined in $(3)$, $q_1$ is the guiding orientation, $\theta$ is the relative angle, and $\langle q_1,q_x \rangle$ is as defined in $(2)$.
I am struggling with how to choose parameters $u_1,u_2,u_3$ to satisfy $(4)$, especially given the need for uniform sampling of the orientation space. The expanded equation $(4)$ is not trivial to simplify. Any help would be greatly appreciated.