This is a quick extension to a previous post on uniform points on disc, circle, sphere and spherical caps^{1}. We can also move away from trig or the rejection method by approximing equalarea disc^{2}. I’ll assume basic understanding of representing rotations by quaternions. For a matrix simply follow through the deriviation by applying the similarity transform to an orthogonal basis.
Informally there are three general types of random rotations in 3D…one, two and three degrees of freedom.
About fixed axis one degree of freedom
The simplest general case is generating a random rotation about a fixed axis (considered as unit bivector $\mathbf{u}$) which is equivalent to the two dimensional case. We can generate a point on the halfcircle (positive in $x$) and set the scalar part of the quaternion to the $x$ coordinate and set the bivector part to $\mathbf{u}$ scaled by the $y$ result.
For nonuniform, if the distribution on the halfcircle has angle probability function $P(\theta)$, then the result probability of the rotation angle is $P(2\theta)$.
With a fixed reference direction two degrees of freedom
In axisangle speak moving to two degress of freedom fixes the axisofrotation to a plane. We can choose a fixed reference direction which is orthogonal to that plane (say positive $\mathbf{z}$), generate a uniform point on the sphere and find the rotation.
First generate a uniform point on the disc: $p=\left(x,~y\right)$ and map to the sphere^{1}:
Find the quaternion that rotates $\left(0,0,1\right)$ to some $\left(x,~y,~z\right)$^{3}:
Composing the previous operations results in:
Breaking this down we have:
 the direction of the axis of rotation: $(y,~x,0)$
 CCW orthogonal to point on disc (which is logically in $\mathbf{xy}$plane)
 magnitude is logically $\sin\left(\frac{\theta}{2}\right)$
 the scalar part is logically $\cos\left(\frac{\theta}{2}\right) = \sqrt{1\sin(\theta)^2} = \sqrt{1d}$
We could have directly generated this result but requires figuring out how to set up a Jacobian…and I hate doing that. Since our values on the disc are random, we can also simply reorder the bivector part and drop the sign (not rotate) without any impact, giving us:
For nonuniform (with the reordering and sign drop) then if the probability on the disc is $P(x,y)$ is the probability of a rotation of $2 \text{ acos}\left(1(x^2+y^2)\right)$ about axis in direction $(x,y,0)$. More generally we can shape the disc distribution and convert to cap.
The Real Deal (arbitrary) three degrees of freedom
George Marsaglia in same paper^{4} as uniform points on the 3D sphere $\left(\mathbb{S}^2\right)$ handwavingly^{5} presented a method for the 4D sphere $\left(\mathbb{S}^3\right)$.
Given two uniform points on the unit disc: $~p_0,~p_1$ and their norm: $d_n=p_n \cdot p_n$:
Making an arbitrary choice for the scalar part ($x_0$) the quaternion is:
References and Footnotes

“Uniform points on disc, circle, sphere and caps”, (local post) ↩ ↩^{2}

“Square/Disc mappings”, (local post) ↩

“Implied normals (unit bivectors as rotations)”, (local post) ↩

“Choosing a Point from the Surface of a Sphere”, George Marsaglia, 1972 (PDF) ↩

Per Vognsen provides a simple justification in this gist: (link) ↩