Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

Direction cone

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.

share|improve this question
    
In what way do you want the random direction vectors to be uniformly distributed? Your bullet points seem to indicate you're trying to make them uniform with respect to what points inside an interior circle they intersect, whereas in my mind a uniform distribution should be created according to the points inside a spherical cap. –  anon Aug 10 '11 at 22:42
    
Yes, and this is why the final step is to normalize the vector, which maps it to a spherical cap. Actually, I'm not even sure if "uniformly distributed" is well-defined in this situation. –  Trillian Aug 10 '11 at 22:49
1  
Uniformity of distribution is not generally preserved after mapping one space to another (e.g. disc to cap). –  anon Aug 10 '11 at 23:36
1  
math.stackexchange.com/questions/44689/… perhaps helps. It is on how to generate a unit vector on a sphere. –  Shiyu Aug 11 '11 at 1:25
2  
@Shiyu Yes, that's by far the simplest solution. Your random direction is obtained by generating a uniform variate in $[0,2\pi)$ for the angle around the axis (azimuth) and another uniform variate in $[\cos(\theta),1]$, which determines the dot product between the random direction and the unit vector along the axis. The two are readily re-expressed as a direction vector. The illustration in your link demonstrates why this is uniformly distributed within the solid angle subtended by the cone. –  whuber Aug 12 '11 at 22:12
add comment

7 Answers

up vote 9 down vote accepted

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.

share|improve this answer
add comment

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

share|improve this answer
    
Thanks for your answer, though honestly I didn't understand much of it. I guess I have some wikipedia reading to do on polar and azimuthal angles. –  Trillian Aug 11 '11 at 0:07
    
The map $\theta=f(x)$ is actually quite simple: Since $z$ is uniformly distributed, $\theta$ can be obtained as $\theta=\arccos z$; see Christian's and my answer. –  joriki Oct 1 '12 at 20:16
add comment

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.

share|improve this answer
add comment

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.

share|improve this answer
    
That looks wise but I can't make much sense of it past the first sentence. I think I should be able to implement it anyways, but I don't see where the original cone axis vector comes into this. –  Trillian Aug 11 '11 at 0:11
    
@Trillian My answer assumed the cone is oriented along $z$ axis in the positive direction. You could apply rotation matrix to the generated vector to orient cone's axes any way you like. I will try to expand the answer to make it more detailed. –  Sasha Aug 11 '11 at 1:48
    
I didn't check whether the $\theta$ sampling is correct, but it's unnecessarily complicated; see Christian's or my answer. –  joriki Oct 1 '12 at 20:14
1  
@joriki While Chrisian's answer is very nice, I do not see mine as unnecessarily complicated. It answers the question posed and certainly not worthy of the downvote somebody awarded it. Christian's answer as well as yours uses the [fact](This question has been posed on math.SE/question/185298) that $z$ component is uniformly distributed. It's truncation to $(1,\cos(\theta))$ yields the desired distribution. What I wrote is exactly the same! My recipe: $$(x,y,z) = (\sin(\vartheta) \cos(\phi), \sin(\vartheta) \sin(\phi), \cos(\vartheta) ) $$ where $\cos(\vartheta) = 1-u + u \cos(\theta)$. –  Sasha Oct 1 '12 at 22:13
    
By as you see $\cos(\vartheta)$ is exactly what you call $z$. I am sorry to see people hastily downvote (as well as upvote) without checking the result for themselves, based on comments like this one. –  Sasha Oct 1 '12 at 22:15
show 3 more comments

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}}$

get a random perpendicular

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.

rotate around that random perpendicular

Result with a 45 degree spread:

45

45

share|improve this answer
2  
This does not give a uniform distribution, because (i) the axis is not uniformly distributed in the plane normal to $\vec v$, and (ii) you need to do something involving $\cos^{-1}$ instead of taking square root; if you let $\theta=\pi$, the square root will oversample the opposite pole. –  Rahul Aug 14 '12 at 23:19
    
Alright, say: i) we choose the random axis to be a uniformly distributed vector on the unit sphere. with regard to ii) I experimented with acos but the distributions weren't correct (interestingly, taking $arccos( \theta^2 )$ takes a $\theta$ sized chunk out of a hemisphere) –  bobobobo Aug 15 '12 at 15:15
add comment

I recommend you should use a rotation matrix for an arbitrary axis in 3D WP has articles : http://en.wikipedia.org/wiki/Rotation_matrix

share|improve this answer
add comment

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.

share|improve this answer
    
This won't result in a uniform distribution, since you're uniformly sampling $\theta$ where you should be uniformly sampling $z$ (see Christian's or my answer). –  joriki Oct 1 '12 at 20:11
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.