Fast uniformly distributed random points on the surface of a unit hemisphere

Question

I am trying to generate uniform random points on the surface of a unit sphere for a Monte Carlo ray tracing program. When I say uniform I mean the points are uniformly distributed with respect to surface area. My current methodology is to calculate uniform random points on a hemisphere pointing in the positive z axis and base in the x-y plane.

The random point on the hemisphere represents the direction of emission of thermal radiation for a diffuse grey emitter.

I achieve the correct result when I use the following calculation :

Note : dsfmt* is will return a random number between 0 and 1.

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt); zenith = asin(sqrt(dsfmt_genrand_close_open(&dsfmtt)));  // Calculate the cartesian point osRay.c._x = sin(zenith)*cos(azimuthal);  osRay.c._y = sin(zenith)*sin(azimuthal); osRay.c._z = cos(zenith); 

However this is quite slow and profiling suggests that it takes up a large proportion of run time. Therefore I sought out some alternative methods:

The Marsaglia 1972 rejection method

do {    x1 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;    x2 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;    S = x1*x1 + x2*x2; } while(S > 1.0f);   osRay.c._x = 2.0*x1*sqrt(1.0-S); osRay.c._y = 2.0*x2*sqrt(1.0-S); osRay.c._z = abs(1.0-2.0*S); 

Analytical cartesian coordinate calculation

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt); u = 2*dsfmt_genrand_close_open(&dsfmtt) -1; w = sqrt(1-u*u);  osRay.c._x = w*cos(azimuthal); osRay.c._y = w*sin(azimuthal); osRay.c._z = abs(u); 

While these last two methods run serval times faster than the first, when I use them I get results which indicate that they are not generating uniform random points on the surface of a sphere but rather are giving a distribution which favours the equator.

Additionally the last two methods give identical final results however I am certain that they are incorrect as I am comparing against an analytical solution.

Every reference I have found indicates that these methods do produce uniform distributions however I do not achieve the correct result.

Is there an error in my implementation or have I missed a fundamental idea in the second and third methods?

Answer 1

The simplest way to generate a uniform distribution on the unit sphere (whatever its dimension is) is to draw independent normal distributions and normalize the resulting vector.

Indeed, for example in dimension 3, e^(-x^2/2) e^(-y^2/2) e^(-z^2/2) = e^(-(x^2 + y^2 + z^2)/2) so the joint distribution is invariant by rotations.

This is fast if you use a fast normal distribution generator (either Ziggurat or Ratio-Of-Uniforms) and a fast normalization routine (google for "fast inverse square root). No transcendental function call is required.

Also, the Marsaglia is not uniform on the half sphere. You'll have more points near the equator since the correspondence point on the 2D disc <-> point on the half sphere is not isometric. The last one seems correct though (however I didn't make the calculation to ensure this).