Generate a random direction within a cone

Question

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.

Answer 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.

Answer 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.

Answer 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

Answer 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.

Answer 5

let's for simplicity generate points within an angle ALPHA of the north pole

P = (1 - cos(ALPHA))/2

U = random(0,1) a random value in the range (0,1)

PHI = acos(1-2*P*U)

THETA = random (-PI,PI)

PHI is the angle from the north pole

Answer 6

So you want to uniformly sample on a spherical cap. With the notations of this link, here is some pseudo-code which performs the sampling on the cap displayed in the above link (then it remains to perform a rotation):

stopifnot(h > 0 && h < 2 * R)
xy = uniform_point_in_unit_circle()
k = h * dotproduct(xy)
s = sqrt(h * (2 * R - k))
x = xy[1]
y = xy[2]
return (s * x, s * y, R - k)

I have this code in my archives, unfortunately I don't remember how I got it.

I put the image of the link here, in case it changes:

Answer 7

I think you have 2 options:

a) Pre-generate a sphere of uniformly distributed random vectors on the sphere described on page 19 of the Total Compendium. When you are "picking" vector from that collection, reject any result with angle greater than the $\theta$ you have indicated in your diagram. You should probably precompute $\cos \theta$ and use a simple dot product to find if the random vector chosen is acceptable or not.

b) Use the following 3 steps to get a random vector that is at most $\theta$ degrees from $\vec{v}$:

1) Get any vector $\vec{\text{random}}$ on the unit sphere.

2) Cross $\vec{\text{random}}$ with $\vec{v}$ to get $\vec{\text{axis}}$

3) Rotate $\vec{v}$ about $\vec{\text{axis}}$ by $\phi$ degrees, where $\phi = \sqrt{\text{random}( 0, \theta^2 )}$ ($\phi$ is the square root of uniformly distributed random scalar between 0 and $\theta^2$). This distribution of $\phi$ is important to avoid crowding at the pole.

Result with a 45 degree spread:

Answer 8

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.