Unit Sphere and Surface Points

[Ganado]

<This post was originally a reply to post that has since been removed. I got some mathematical details wrong in my original reply, see Thomas' and lastchance's replies>

A unit sphere has a radius 1, and you can break a point on the surface on the sphere down as a vector using spherical coordinates,
http://tutorial.math.lamar.edu/Classes/CalcIII/SphericalCoords.aspx

x = r sin(phi) cos(theta)
y = r sin(phi) sin(theta)
z = r cos(phi)

Now all you need to do is uniformly generate a random phi between 0 and pi, and a random theta between 0 and 2 pi. You can generate a random floating-point number in modern C++ using the <random> library with uniform_real_distribution.
EDIT: Look at posts below. It may be more complicated than just uniform phi/theta.

https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution
http://www.cplusplus.com/reference/random/uniform_real_distribution/operator()/

#include <random>
#include <iostream>
#include <chrono>
 
int main()
{
    const double Pi = 3.14159265359; // ...and the rest
    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
    std::mt19937 gen(seed);
    std::uniform_real_distribution<double> theta_dist(0.0, 2.0 * Pi);
    
    for (int n = 0; n < 10; ++n) {
        std::cout << theta_dist(gen) << '\n';
    }
    std::cout << '\n';
}

This is seeded based on the current system clock (in seconds since epoch time).

[Thomas Huxhorn]

I have this fortran code to pick a random point in the unit ball in 3D. With some explanation.

https://people.sc.fsu.edu/~jburkardt/f_src/geometry/geometry.f90

!*****************************************************************************80
!
!! BALL01_SAMPLE_3D picks a random point in the unit ball in 3D.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    26 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) SEED, a seed for the random 
!    number generator.
!
!    Output, real ( kind = 8 ) P(3), the sample point.
!
  implicit none

  integer ( kind = 4 ), parameter :: dim_num = 3

  real ( kind = 8 ) p(dim_num)
  real ( kind = 8 ) phi
  real ( kind = 8 ) r
  real ( kind = 8 ) r8_acos
  real ( kind = 8 ), parameter :: r8_pi = 3.141592653589793D+00
  integer ( kind = 4 ) seed
  real ( kind = 8 ) theta
  real ( kind = 8 ) u(dim_num)
  real ( kind = 8 ) vdot

  call r8vec_uniform_01 ( dim_num, seed, u )
!
!  Pick a uniformly random VDOT, which must be between -1 and 1.
!  This represents the dot product of the random vector with the Z unit vector.
!
!  Note: this works because the surface area of the sphere between
!  Z and Z + dZ is independent of Z.  So choosing Z uniformly chooses
!  a patch of area uniformly.
!
  vdot = 2.0D+00 * u(1) - 1.0D+00

  phi = r8_acos ( vdot )
!
!  Pick a uniformly random rotation between 0 and 2 Pi around the
!  axis of the Z vector.
!
  theta = 2.0D+00 * r8_pi * u(2)
!
!  Pick a random radius R.
!
  r = u(3) ** ( 1.0D+00 / 3.0D+00 )

  p(1) = r * cos ( theta ) * sin ( phi )
  p(2) = r * sin ( theta ) * sin ( phi )
  p(3) = r * cos ( phi )

  return
end

The routine BALL01_SAMPLE_2D is also very interesting.

[Ganado]

Nice link, but I'm confused why the dot product (and acos calculation) is needed in your example?
It's very possible I'm missing something essential here.
Edit: Oh I think the difference might be that example is uniformly spreading out the volume.

picks a random point in the unit ball in 3D.

[lastchance]

Is Z the coordinate along the vertical axis? There is more surface area of the sphere near the equator (z=0) than near the poles (z=1 or -1). Surely you need to weight your probabilities in inverse proportion to the length of a line of latitude.

[Ganado]

I agree, there is unequal surface area in a horizontal slice near the equator vs. near the top, but z isn't being generated uniformly, theta/phi is.

I think it's as simple as uniformly generating theta/phi and then applying spherical coordinates... (if we're talking about the surface of a sphere, not inside of it) but show me if I'm wrong, I could be wrong.
[If anyone wanted to do it, this could be checked numerically just by looking at the density of generated points at each point on the surface]

[lastchance]

It isn't uniform in colatitude and longitude. (Sorry, I have to use those terms because I use spherical coordinates with theta and phi the opposite way round to you. Colatitude is angle from the north pole. Longitude is angle from one meridion.)
A small increment of surface on the sphere r=1 is
sin(colatitude) d(colatitude) d(longitude)
or
d(-cos(colatitude)) d(longitude)
or
d(-z) d(longitude)

Ah, I see what @Thomas Huxhorn means now - it is uniform in z and longitude, with z ranging from -1 to 1 and longitude from 0 to 2.pi, so giving a total surface area of 4 pi (which is correct).

So, generate z uniform in [-1,1]
and longitude uniform in [0, 2 pi]

You can recover colatitude from cos(colatitude)=z, which is where the acos() function is coming from.


Grrr. I code in Fortran most of the time but it took me too long to decipher that. Note that Fortran's standard random-number generation is uniform on [0,1) (although it's useful that you can get an entire arrayful at once, rather than repeated calls).



Here's an attempt in C++.

#include <iostream>
#include <iomanip>
#include <random>
#include <chrono>
#include <cmath>
using namespace std;


const double PI = 3.141592653589793238;
unsigned seed = chrono::system_clock::now().time_since_epoch().count();
mt19937 gen( seed );
uniform_real_distribution<double> dist( 0.0, 1.0 );


//====================================================================
// Return random point on the unit sphere r=1
// Note that it returns both the cartesian coordinates x, y, z
// and the spherical angles colatitude and longitude (in radians)
//====================================================================
void getPoint( double &x, double &y, double &z, double &colatitude, double &longitude )
{
   z = -1.0 + 2.0 * dist( gen );       // uniformly distributed z
   longitude = 2.0 * PI * dist( gen ); // uniformly distributed longitude

   colatitude = acos( z );
   double rh = sqrt( 1.0 - z * z );    // horizontal radius (= r sin(colatitude), with r=1 )
   x = rh * cos( longitude );
   y = rh * sin( longitude );
}


int main()
{
   #define SP << setw(10) << setprecision(4) << fixed <<
   const int N = 10;                        
   double x, y, z, colatitude, longitude;
   
   cout SP "x" SP "y" SP "z" << '\n';
   for ( int i = 0; i < N; i++ )
   {
      getPoint( x, y, z, colatitude, longitude );
      cout SP x SP y SP z << '\n';
   }
}

         x         y         z
   -0.8246    0.1450   -0.5469
   -0.2953    0.0332   -0.9548
    0.9162   -0.0761    0.3934
    0.2246    0.5299   -0.8178
   -0.6054   -0.2804    0.7449
   -0.5632    0.6979    0.4424
   -0.3306   -0.7446   -0.5799
    0.3448   -0.7311   -0.5888
   -0.6028   -0.0653    0.7952
   -0.0447    0.6338    0.7722

[Ganado]

Thanks for the reply, I will have to study this.

And I should have been more careful about the phi/theta terminology, I know some textbooks/fields use the opposite notation... ugh...

[Duthomhas]

Erg, anything using spherical coordinates is biased toward the poles, because the two randomly-generated numbers are not independent.

Generate a 3D point in the unit square (3 RNGs, which are independent), toss all points outside the sphere (otherwise you will again introduce bias towards the corners), then simply normalize the vector to unit length.

To avoid some FP issues you may also want to toss all points creating a vector less than .01 unit length (or some reasonable fraction) before normalizing.

Hope this helps.

[lastchance]

Yes, that would work, but then you are generating 3 random numbers for each point (and then rejecting a sizeable fraction outside the sphere), whereas @Thomas Huxhorn's method (reduced to 2d, because we are looking for surface rather than volume, fixing r=1) would only need 2 rng calls and all would be used.

It is uniform and independent in the vertical coordinate (z) and longitudinal angle.

Wish there was consensus about the order of phi and theta, though! My preferred notation (opposite to the earlier part of the thread) as here:
https://math.stackexchange.com/questions/131735/surface-element-in-spherical-coordinates
with the formulation in terms of z and longitude in
https://math.stackexchange.com/questions/131735/surface-element-in-spherical-coordinates/870413#870413

See also:
http://mathworld.wolfram.com/SpherePointPicking.html



EDIT: LOOKS LIKE THE ORIGINAL POST BY @GilbertJobs MAY HAVE BEEN TROLLING (interesting question though it was) AND IS PROBABLY COPIED FROM REDDIT.

[coder777]

Look at the posted links. OP is most likely spam.

[Ganado]

Gotta be kidding me, it edited the post. Can anyone remove its post?

That wolfram article explains it pretty well, I see why now.

[lastchance]

I've reported (and thence removed) the OP's post with the spam links that he/she edited in. However, it was quite an educative thread ... which you are now proud "owner" of, @Ganado!