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.

Solution 1:

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.

Solution 2:

Your first four bulleted points are absolutely correct as stated. So you now have two mutually orthogonal unit vectors ${\bf u}$, ${\bf v}$, both of them orthogonal to the given axis ${\bf a}$ of the cone, where $|{\bf a}|=1$. The random unit vector within the cone will be a vector ${\bf x}$ of the form $${\bf x}=\sin\theta\bigl(\cos\phi\,{\bf u}+\sin\phi\, {\bf v}\bigr)+\cos\theta\,{\bf a}\ .$$ Here $\phi$ has to be chosen uniformly in $[-\pi,\pi]$, and $\theta$ has to be chosen according to some distribution yet to be determined in the interval $[0,\theta_0]$, where $\theta_0$ is the angle denoted by $\theta$ in your figure. Note that $\theta_0$ is restricted to the interval $[0,\pi]$, not to $[0,2\pi]$ as indicated in your question.

Now we want the vectors ${\bf x}$ to be equidistributed on the spherical cap given by $0\leq\theta\leq\theta_0$, where "equidistribution" refers to the area measure on the sphere $S^2$.

Here the following elementary fact comes to our help: When a cylinder $C$ of height $2$ is wrapped around the equator of a unit sphere then any two planes orthogonal to the axis of the cylinder determine a "spherical annulus" and a cylindrical annulus, and the areas of these two annuli coincide. This implies that points uniformly distributed on the cylinder $C$ "project" to points uniformly distributed on the sphere $S^2$.

This observation boils down to the following recipe: Choose $z$ uniformly distributed on the interval $[\cos\theta_0,1]$, and put $\theta:=\arccos(z)$. Furthermore let $\phi$ be uniformly distributed on $[0,2\pi]$. Then the point $${\bf x}=\sin\theta\bigl(\cos\phi\,{\bf u}+\sin\phi\, {\bf v}\bigr)+\cos\theta\,{\bf a}$$ will be uniformly distributed on the spherical cap in question.

Solution 3:

I'm assuming that the kind of uniform distribution you want is one in which probability is proportional to solid angle.

I would generate a polar angle $\theta$ and an azimuthal angle $\phi$, and then determine the vector from that. $\phi$ can be uniformly distributed from 0 to $2\pi$.

You do not want $\theta$ to be uniformly distributed, however; you want its probability distribution to be proportional to $\sin\theta$, since that is how the solid angle varies with $\theta$.

One way to accomplish this is to generate a value for $\theta$ with a uniform distribution between 0 and $\theta_{max}$, but then accept it only with probability $\sin\theta$; if it's not accepted, you try again. This has the disadvantage that it can be inefficient, especially if $\theta_{max}$ is small.

A better technique is to find a map $\theta=f(x)$ such that if $x$ is uniformly distributed on [0,1], $\theta$ is distributed on $[0,\theta_{max}]$ with a probability distribution proportional to $\sin\theta$. This can be done by using the general relationship between the probability distribution of two random variables that are related by some function; it's essentially just an application of the chain rule. This is a standard technique.

WP has articles on both techniques:

http://en.wikipedia.org/wiki/Rejection_sampling

http://en.wikipedia.org/wiki/Inverse_transform_sampling

Solution 4:

Assume the direction of cone's axis is given by $(0,0,1)$.

The random direction within the code corresponds to a random point on a sphere, truncated down to the part cut out by a cone. The measure on a sphere in spherical coordinates factors $\cos \vartheta d \theta d \varphi$. The cone corresponds to a region given by $ 0\leqslant \vartheta\leqslant\theta$ and $0\leqslant \varphi <2 \pi$.

Hence we should draw $0 \leqslant \vartheta < \theta$ from density $$f(\vartheta) \,\mathrm{d}\vartheta = \frac{1}{1-\cos \theta} \sin \vartheta \, \mathrm{d} \vartheta = \mathrm{d} \left( \frac{1-\cos \vartheta }{1-\cos \theta}\right) = \mathrm{d} F(\vartheta)$$ then uniform $\varphi$ from $[0, 2\pi)$, and then build $$(x,y,z) =(\sin \vartheta \cos\varphi, \sin \vartheta \sin\varphi, \cos \vartheta)$$

To sample $\vartheta$ we compute quantile of the distribution for $\vartheta$ and apply it to the uniform variate (see inversion method): $$ \vartheta = F^{(-1)}(u) = \arccos( 1-u + u \cos \theta) $$ where $u$ follows uniform distribution on unit interval. Notice that $z$-component $$ \cos \vartheta = 1 \cdot (1-u) + \cos\theta \cdot u $$ is uniformly distributed on $(1,\cos\theta)$, in agreement with the algorithm of Christian and Joriki. Also see this relevant question.