Generate a random direction within a cone I have a normalized $3D$ vector giving a direction and an angle that forms a cone around it, something like this:



I'd like to generate a random, uniformly distributed normalized vector for a direction within that cone. I would also like to support angles greater than pi (but lower or equal to $2\pi$), at which point the shape becomes more like a sphere from which a cone was removed. How can I proceed?

I thought about the following steps, but my implementation did not seem to work:


Find a vector normal to the cone axis vector (by crossing the cone axis vector with the cardinal axis that corresponds with the cone axis vector component nearest to zero, ex: $[1 0 0]$ for $[-1 5 -10]$)
Find a second normal vector using a cross product
Generate a random angle between $[-\pi, \pi]$
Rotate use the two normal vectors as a $2D$ coordinate system to create a new vector at the angle previously generated
Generate a random displacement value between $[0, \tan(\theta)]$ and square root it (to normalize distribution like for points in a circle)
Normalize the sum of the cone axis vector with the random normal vector times the displacement value to get the final direction vector


[edit] After further thinking, I'm not sure that method would work with theta angles greater or equal to pi. Alternative methods are very much welcome.
I'm surprised how many bad, suboptimal and/or overcomplicated answers this question has inspired, when there's a fairly simple solution; and that the only answer that mentions and uses the most relevant fact, Christian Blatter's, didn't have a single upvote before I just upvoted it.

The $2$-sphere is unique in that slices of equal height have equal surface area. That is, to sample points on the unit sphere uniformly, you can sample $z$ uniformly on $[-1,1]$ and $\phi$ uniformly on $[0,2\pi)$. If your cone were centred around the north pole, the angle $\theta$ would define a minimal $z$ coordinate $\cos\theta$, and you could sample $z$ uniformly on $[\cos\theta,1]$ and $\phi$ uniformly on $[0,2\pi)$ to obtain the vector $(\sqrt{1-z^2}\cos\phi,\sqrt{1-z^2}\sin\phi,z)$ uniformly distributed as required.

So all you have to do is generate such a vector and then rotate the north pole to the centre of your cone. If the cone is already thus centred, you're done; if it's centred on the south pole, just invert the vector; otherwise, take the cross product of the cone's axis with $(0,0,1)$ to get the direction of the rotation axis, and the scalar product to get the cosine of the rotation angle. Or if you prefer you can apply your idea of generating two orthogonal vectors, in the manner Christian described.
I've thought of another approach that seems to produce good results from quick visualization tests. I'd love to get feedback on this.

The 2D case is much simpler, I just have to add an angle within [-theta, theta] to the "cone" axis vector's angle and create a vector from the new angle. The 3D case is like a solid of revolution. I can generate an angle within [0, theta] and then move it in 3D space to a random rotation about the cone axis vector. I end up with the following algorithm:

axis = cone axis vector
angleAroundAxis = random(0, 2 * pi)
angleAwayFromAxis = random(0, theta)

matrix1 = rotation matrix of 'angleAroundAxis' around 'axis'
anyNormal = any vector that is normal to 'axis'

normal = anyNormal * matrix1
binormal = normal x axis

matrix2 = rotation matrix from 'axis', 'normal' and 'binormal' axes

result = [cos(angleAwayFromAxis) sin(angleAwayFromAxis) 0] * matrix2


Sorry if this is using non-standard notation, I'm not very familiar with formal math notation.