Computing Random Points on Sphere
Dear, I'm a newcomer on the boost list. I subscribed to share the results of my experiments on picking random points on sphere. This procedure is often a bottleneck for very fast Monte-Carlo simulations of physical process where you must pick vectors in random directions. I compiled a document that explain the experiments and the results. The conclusion contains some propositions to modify boost::random_on_sphere and I would like to know if a new implementation following the presented ideas would be welcome. Colas Schretter
I doubt that the trig-version gives a good distribution: inlineVectortrig(uniform_01<generator_type>&random) { constfloatphi=2*M_PI*random(); constfloatcos_theta=2*random()-1; constfloatsin_theta=sqrt(1-cos_theta*cos_theta); returnVector(cos_theta,sin_theta*cos(phi),sin_theta *sin(phi)); } The result ought to be more dense at the poles. /$ 2007/6/5, more effective thinking in the exceptional C++ programming language <effective.thinking@gmail.com>:
Dear,
I'm a newcomer on the boost list. I subscribed to share the results of my experiments on picking random points on sphere. This procedure is often a bottleneck for very fast Monte-Carlo simulations of physical process where you must pick vectors in random directions. I compiled a document that explain the experiments and the results. The conclusion contains some propositions to modify boost::random_on_sphere and I would like to know if a new implementation following the presented ideas would be welcome.
Colas Schretter
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
Henrik Sundberg wrote:
I doubt that the trig-version gives a good distribution:
inlineVectortrig(uniform_01<generator_type>&random) { constfloatphi=2*M_PI*random(); constfloatcos_theta=2*random()-1; constfloatsin_theta=sqrt(1-cos_theta*cos_theta); returnVector(cos_theta,sin_theta*cos(phi),sin_theta *sin(phi)); }
The result ought to be more dense at the poles.
The posted trig method is correct, it chooses a horizontal cutting plane uniformly. You can show that the surface area of a sphere contained between two z-coordinates only depends on the difference between the two coordinates, not their location, so this generates the correct distribution. However, the "naive" method is incorrect, the resulting points are not uniformly distributed over the sphere. The correct version of the "naive" method is given here: http://mathworld.wolfram.com/SpherePointPicking.html -Lewis
The posted trig method is correct, it chooses a horizontal cutting plane uniformly. You can show that the surface area of a sphere contained between two z-coordinates only depends on the difference between the two coordinates, not their location, so this generates the correct distribution.
However, the "naive" method is incorrect, the resulting points are not uniformly distributed over the sphere. The correct version of the "naive" method is given here: http://mathworld.wolfram.com/SpherePointPicking.html
I take it back... this method does work OK after all. I believe those listed at the above link are superior, though. -Lewis
Lewis Hyatt wrote:
Henrik Sundberg wrote:
I doubt that the trig-version gives a good distribution:
inlineVectortrig(uniform_01<generator_type>&random) { constfloatphi=2*M_PI*random(); constfloatcos_theta=2*random()-1; constfloatsin_theta=sqrt(1-cos_theta*cos_theta); returnVector(cos_theta,sin_theta*cos(phi),sin_theta *sin(phi)); }
The result ought to be more dense at the poles.
The posted trig method is correct, it chooses a horizontal cutting plane uniformly. You can show that the surface area of a sphere contained between two z-coordinates only depends on the difference between the two coordinates, not their location, so this generates the correct distribution.
However, the "naive" method is incorrect, the resulting points are not uniformly distributed over the sphere. The correct version of the "naive" method is given here: http://mathworld.wolfram.com/SpherePointPicking.html
-Lewis
_______________________________________________
Also see here for a discussion for 4 different techniques (and a proof for the trig version!): http://www.math.niu.edu/~rusin/known-math/96/sph.rand Jeff Faust
2007/6/5, Lewis Hyatt <lhyatt@princeton.edu>:
Henrik Sundberg wrote:
The result ought to be more dense at the poles.
The posted trig method is correct, it chooses a horizontal cutting plane uniformly. You can show that the surface area of a sphere contained between two z-coordinates only depends on the difference between the two coordinates, not their location, so this generates the correct distribution.
Aah, sorry for the noise! /$
Thank you for the involved discussions. Picking random points on the surface of a sphere seems to be a never-ending topic of discussion since ages! Thank to the very interesting post of Jeffrey, I discovered a fourth technique described in the post of George Marsaglia pasted bellow (technique 2). It can be viewed as a combination of rejection sampling and trigonometric relations. It appears to be a clear winner in the point of view of performances. So please consider the following updated document in attachment as applicable. At the end of his post, George Marsalia announce using a very fast Gaussian distributed generator of random numbers. This algorithm could be the key to improve the performances of both the random_on_sphere and the Gaussian distributions of boost. Colas Schretter From: George Marsaglia <geo@stat.fsu.edu> Subject: Random points on a sphere. Date: Tue, 15 Jun 1999 18:03:07 -0400 Newsgroups: sci.math.num-analysis,sci.math,sci.stat.math Questions on methods for random sampling from the surface of a sphere keep coming up, and seemingly get lost as new generations of computer users encounter the problem. I described the methods that I had used for years in "Sampling from the surface of a sphere", Ann Math Stat, v43, 1972, recommending the two that seemed the fastest: 1) Choose standard normal variates x1,x2,x3, put R=1/sqrt(x1^2+x2^2+x3^2) and return the point (x1*R,x2*R,x3*R). 2) Generate u and v, uniform in [-1,1] until S=u^2+v^2 < 1, then return (1-2*S,u*r,v*r), with r=2*sqrt(1-S). The second method is fast and easiest to program---unless one already has a fast normal generator in his library. Method 1) seems by far the fastest (and easily applies to higher dimensions) if my latest normal generator is used. It produces normal variates at the rate of seven million per second in a 300MHz Pc. George Marsaglia On 6/5/07, Henrik Sundberg <storangen@gmail.com> wrote:
2007/6/5, Lewis Hyatt <lhyatt@princeton.edu>:
Henrik Sundberg wrote:
The result ought to be more dense at the poles.
The posted trig method is correct, it chooses a horizontal cutting plane uniformly. You can show that the surface area of a sphere contained between two z-coordinates only depends on the difference between the two coordinates, not their location, so this generates the correct distribution.
Aah, sorry for the noise!
/$ _______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
participants (4)
-
Henrik Sundberg -
Jeffrey Faust -
Lewis Hyatt -
more effective thinking in the exceptional C++ programming language