It's not clear to me what you mean by restricting the rotations to a range $[-\theta,\theta]$. The sign of the rotation angle of a rotation isn't meaningful, since you can instead view the rotation as a rotation about the opposite axis through the negative angle. Thus I'll assume that you intended the rotation angles to be restricted to $[0,\theta]$.
Also, I'll assume that you want the distribution to be uniform with respect to the Haar measure of $SO(3)$, which is what the method you quote produces. You're right that this isn't achieved by a uniform distribution over the angles. The probability of an angle should be proportional to the corresponding slice of the surface of $S^3$, the three-dimensional hyper-sphere in four dimensions parameterizing the unit quaternions, and this is proportional to $\sin^2\alpha\,\mathrm d\alpha$.
You can use the method you describe and then reject rotations with angles you don't want. I suspect that this would be one of the more efficient methods if your angle range is large. As your cutoff $\theta$ decreases, however, this becomes increasingly wasteful, since small angles are unlikely, but in this case it becomes increasingly efficient to use rejection sampling with the $\sin^2\alpha$ distribution approximated by $\alpha^2$.
Thus, you could generate a sample from the density $\alpha^2$ using $u=\alpha^3$ and $\mathrm du=3\alpha^2\mathrm d\alpha$ by generating a uniform sample $u\in[0,\theta^3]$, computing $\alpha=\sqrt[3]u$ and accepting the sample with probability $\sin^2\alpha/\alpha^2$, which will be close to $1$ for small $\theta$.
[Edit in response to the comment:]
Rejection sampling for a density $f(\alpha)$ works by sampling from a different density $g(\alpha)$ and then accepting the samples with probability $f(\alpha)/Mg(\alpha)$, with $M=\sup_\alpha f(\alpha)/g(\alpha)$. In the present case, $f(x)\propto\sin^2\alpha$, $g(x)\propto\alpha^2$ and $\sup_\alpha\sin^2\alpha/\alpha^2=1$, so you can accept with probability $\sin^2\alpha/\alpha^2$, which is greater than $99.7\%$ for angles below $5^\circ$, that is, you almost never have to reject a sample and the sampling is very efficient.
To generate samples with density $g(\alpha)\propto\alpha^2$, you can transform to $u=\alpha^3$, since then we have $\mathrm du=3\alpha^2\mathrm d\alpha$, so $\alpha$ has density $\propto\alpha^2$ if $u$ has constant density, i.e. is uniformly distributed. Thus, you can generate uniform samples on $[0,\theta^3]$ (which would be $u=\theta\cdot\theta\cdot\theta\cdot\operatorname{rand}(1.0)$ in your code), then calculate $\alpha=\sqrt[3]u$, then generate another variable $t$ uniformly distributed on $[0,1]$ and accept $\alpha$ if $t\lt\sin^2\alpha/\alpha^2$. If $\alpha$ is rejected, you try again.
This gives you a rotation angle $\alpha$ with the right distribution; then you just have to generate a random axis and form the corresponding rotation.