Uniform random points on sphere

From CGAFaq

(Redirected from Random Points On Sphere)

There are several methods. Perhaps the easiest to understand is the rejection method: generate random points in an origin- centered cube with opposite corners (r,r,r) and (−r,−r,−r). Reject any point p that falls outside of the sphere of radius r. Scale the vector to lie on the surface. Because the cube to sphere volume ratio is ^π⁄₆, the average number of iterations before an acceptable vector is found is ⁶⁄_π = 1.90986. This essentially doubles the effort, and makes this method slower than the trig method. A timing comparison conducted by Ken Sloan showed that the trig method runs in about ²⁄₃’s the time of the rejection method. He found that methods based on the use of normal distributions are twice as slow as the rejection method.

The trig method. This method works only in 3-space, but it is very fast. It depends on the slightly counterintuitive fact (see proof below) that each of the three coordinates is uniformly distributed on [−1,1] (but the three are not independent, obviously). Therefore, it suffices to choose one axis (Z, say) and generate a uniformly distributed value on that axis. This constrains the chosen point to lie on a circle parallel to the X-Y plane, and the obvious trig method may be used to obtain the remaining coordinates.

Choose z uniformly distributed in [−1,1].
Choose t uniformly distributed on [0, 2 π).
Let r = √(1−z²).
Let x = r * cos(t).
Let y = r * sin(t).

This method uses uniform deviates (faster to generate than normal deviates), and no set of coordinates is ever rejected.

Here is a proof of correctness for the fact that the z-coordinate is uniformly distributed. The proof also works for the x- and y- coordinates, but note that this works only in 3-space.

The area of a surface of revolution in 3-space is given by

A = 2 * π * ∫_a^b f(x) * √(1 + [f ′(x)]²) dx

Consider a zone of a sphere of radius R. Since we are integrating in the z direction, we have

f(z) = √(R² − z²)

f ′(z) = −z / √(R²−z²)

1 + [f ′(z)]² = r² / (R²−z²)

A	=	2 * π * ∫_a^b √(R²−z²) * R/√(R²−z²) dz
	=	2 * π * R ∫_a^b dz
	=	2 * π * R * (b−a)
	=	2 * π * R * h,

where h = b−a is the vertical height of the zone. Notice how the integrand reduces to a constant. The density is therefore uniform.

Here is C code implementing the trig method:

void SpherePoints(int n, double X[], double Y[], double Z[])
{
int i;
double x, y, z, w, t;

for( i=0; i< n; i++ ) {
  z = 2.0 * drand48() - 1.0;
  t = 2.0 * M_PI * drand48();
  w = sqrt( 1 - z*z );
  x = w * cos( t );
  y = w * sin( t );
  printf("i=%d:  x,y,z=\t%10.5lf\t%10.5lf\t%10.5lf\n", i, x,y,z);
  X[i] = x; Y[i] = y; Z[i] = z;
  }
}

A complete package is available at

ftp://cs.smith.edu/pub/code/sphere.tar.gz (4K),

reachable from http://cs.smith.edu/~orourke/ .

If one wants to generate the random points in terms of longitude and latitude in degrees, these equations suffice:

Longitude = random * 360° − 180°

Latitude = [arccos( random * 2 − 1 )/π ] * 180° − 90°

[edit]