22

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?

cubiclewar
  • 1,569
  • 3
  • 20
  • 37
  • 1
    What kind of precision do you need? It might be fastest to use your first algorithm and just explore ways to speed up the trigonometric functions. For example, you could precompute a table of sines and cosines and then use some kind of fast lookup and interpolation instead of a full-precision calculation. I'm not sure how the roundoff error would affect your results, however. – Daniel Pryden Sep 02 '11 at 07:18
  • I have thought about using fast trig functions however I fear the simulation is fairly angle sensitive. Precomputing a table of trig values would also lend to an increase in cache miss. I starting to think the first method might not be uniform and that is the difference. – cubiclewar Sep 02 '11 at 07:29
  • 1
    Did you try to optimize the first method? calling `sin(asin(x))` is equivalent to using `x` with possible sign changes. When you know `cos(z)`, you can calculate `sin(z)` from the equality cos(z)^2 + sin(z)^2 = 1. – quant_dev Sep 02 '11 at 07:46
  • The first method can be slightly optimised however it still will not come close to the speed of the second and third methods. – cubiclewar Sep 02 '11 at 08:20
  • Your first method is not uniform. See my answer below. – TonyK Sep 02 '11 at 08:39
  • To use Marsaglia correctly, you need to generate three uniform coordinates, not two (i.e. a random point in the unit cube), and accept the point only if it lies within the unit sphere. – TonyK Sep 02 '11 at 08:41
  • TonyK I have been thinking that the first method was not uniform, however there it is commonly referenced for this problem that you need a uniform distribution so I assumed it was uniform (my mistake). Any chance you know what type of point distribution it will give? – cubiclewar Sep 02 '11 at 08:43
  • In your first method, `z` is distributed as `cos(asin(sqrt(h)))`, where `h` is uniformly distributed on `[0,1]`. This simplifies to `sqrt(1-h)`, which will skew the distribution heavily away from the equator. Where did you get that formula? – TonyK Sep 02 '11 at 09:13
  • It was presented in a paper circa 1968 called "Application of Monte Carlo to Heat Transfer Problems" by John Howell. I believe it is derived from the total hemispherical emissivity and simplified for a diffuse gray emitter. From the results it definitely skew's the results away from the equator. I think I will have to look into why this is the case. – cubiclewar Sep 02 '11 at 12:04
  • https://github.com/ampl/gsl/blob/644e768630841bd085cb7121085a688c4ff424d0/randist/sphere.c#L66 – zwcloud Jan 13 '19 at 07:40

7 Answers7

16

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

Alexandre C.
  • 55,948
  • 11
  • 128
  • 197
  • Definitely agree. This is my go-to method for uniform points on a sphere. It also generalises to any number of dimensions, if you ever care about that... – Ian Ross Sep 02 '11 at 08:13
  • Indeed. The generalization to arbitrary n is useful, and since when n gets large funny things happen, you have no choice. – Alexandre C. Sep 02 '11 at 08:15
  • Ok I might also consider this, I did come across it in my research however as I will never look at anything above a 2-sphere I dismissed it. – cubiclewar Sep 02 '11 at 08:23
  • @TDawg: you can simply flip the z component if it has wrong sign. – Alexandre C. Sep 02 '11 at 09:06
4

If you take a horizontal slice of the unit sphere, of height h, its surface area is just 2 pi h. (This is how Archimedes calculated the surface area of a sphere.) So the z-coordinate is uniformly distributed in [0,1]:

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
osRay.c._z = dsfmt_genrand_close_open(&dsfmtt);

xyproj = sqrt(1 - osRay.c._z*osRay.c._z);
osRay.c._x = xyproj*cos(azimuthal); 
osRay.c._y = xyproj*sin(azimuthal);

Also you might be able to save some time by calculating cos(azimuthal) and sin(azimuthal) together -- see this stackoverflow question for discussion.

Edited to add: OK, I see now that this is just a slight tweak of your third method. But it cuts out a step.

Community
  • 1
  • 1
TonyK
  • 16,761
  • 4
  • 37
  • 72
3

This should be quick if you have a fast RNG:

// RNG::draw() returns a uniformly distributed number between -1 and 1.

void drawSphereSurface(RNG& rng, double& x1, double& x2, double& x3)
{
    while (true) {
        x1 = rng.draw();
        x2 = rng.draw();
        x3 = rng.draw();
        const double radius = sqrt(x1*x1 + x2*x2 + x3*x3);
        if (radius > 0 && radius < 1) {
            x1 /= radius;
            x2 /= radius;
            x3 /= radius;
            return;
        }
    }   
}

To speed it up, you can move the sqrt call inside the if block.

quant_dev
  • 6,181
  • 1
  • 34
  • 57
  • 1
    This is not uniform on the sphere. It gives more mass to the point near the diagonal (regions corresponding to the vertices of the unit cube). – Alexandre C. Sep 02 '11 at 08:07
  • 3
    Did you notice the check for the condition `radius < 1`? – quant_dev Sep 02 '11 at 08:11
  • 1
    This seems like a variant of the Marsaglia method but runs a fair bit slower and gives the same result. – cubiclewar Sep 02 '11 at 09:27
  • interesting. So with the radius > 0 and radius < 1 we are truncating the cube point cloud to be a sphere point cloud. Geinius that works! – cmarangu Apr 06 '21 at 19:39
2

Have you tried getting rid of asin?

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
sin2_zenith = dsfmt_genrand_close_open(&dsfmtt);
sin_zenith = sqrt(sin2_zenith);

// Calculate the cartesian point
osRay.c._x = sin_zenith*cos(azimuthal); 
osRay.c._y = sin_zenith*sin(azimuthal);
osRay.c._z = sqrt(1 - sin2_zenith);
Pablo
  • 8,644
  • 2
  • 39
  • 29
  • This gives the correct result as you would expect (it is a manipulation using trigonometric identifies) and give a decent run time reduction. However it is still not in the same league as the second method. – cubiclewar Sep 02 '11 at 08:38
1

I think the problem you are having with non-uniform results is because in polar coordinates, a random point on the circle is not uniformly distributed on the radial axis. If you look at the area on [theta, theta+dtheta]x[r,r+dr], for fixed theta and dtheta, the area will be different of different values of r. Intuitivly, there is "more area" further out from the center. Thus, you need to scale your random radius to account for this. I haven't got the proof lying around, but the scaling is r=R*sqrt(rand), with R being the radius of the circle and rand begin the random number.

carlpett
  • 12,203
  • 5
  • 48
  • 82
  • The OP is trying to simulate a uniform distribution on the surface of a 3D sphere, not on the area of a 2D circle. – quant_dev Sep 02 '11 at 07:44
  • @quant_dev: Yes, but as I read his code, he is doing it by using polar coordinates and then calculating a height above the plane. So unless I'm reading it wrong, my answer is still relevant – carlpett Sep 02 '11 at 07:47
  • I understand that the area fractions change with the zenith angle which is why the zenith angle cannot be randomly selected in the range (0,PI/2) for a hemisphere but I am not sure of the implications on the proposed functions. I thought they already accounted for the changing surface area. – cubiclewar Sep 02 '11 at 07:54
1

The second and third methods do in fact produce uniformly distributed random points on the surface of a sphere with the second method (Marsaglia 1972) producing the fastest run times at around twice the speed on an Intel Xeon 2.8 GHz Quad-Core.

As noted by Alexandre C there is an additional method using the normal distribution which expands to n-spheres better than the methods I have presented.

This link will give you further information on selecting uniformly distributed random points on the surface of a sphere.

My initial method as pointed out by TonyK does not produce uniformly distributed points and rather bias's the poles when generating the random points. This is required by the problem I am trying to solve however I simply assumed it would generate uniformly random points. As suggested by Pablo this method can be optimised by removing the asin() call to reduce run time by around 20%.

cubiclewar
  • 1,569
  • 3
  • 20
  • 37
  • Did you test the method with the unit normal distributed variables ? I believe that it should be competitive with Marsaglia's method provided you use a fast normal generator (ratio of uniforms or Ziggurat). – Alexandre C. Sep 05 '11 at 14:32
  • Alexandre, no I haven't tried it because I found that due to the property I am trying to model I actual need to use a distribution which is not uniform (instead concentrated at the poles). However I do have another area where I could apply this and will test it soon. – cubiclewar Sep 06 '11 at 11:55
-1

1st try (wrong)

point=[rand(-1,1),rand(-1,1),rand(-1,1)];
len = length_of_vector(point);

EDITED:

What about? 

while(1)
 point=[rand(-1,1),rand(-1,1),rand(-1,1)];
 len = length_of_vector(point);
 if( len > 1 )
     continue;
 point = point / len
     break

Acception is here approx 0.4. Than mean that you will reject 60% of solutions.

Luka Rahne
  • 10,336
  • 3
  • 34
  • 56