Generating uniform random numbers on a sphere/circle filled with a cube/square

Question

I am trying to generate random points on a sphere that is filled with a cube. Because I had no idea how to do that i started with 2d. (A circle filled with a quadrat.)

What I am trying to do: Generating random points inside the outer circle, but outside the green square.

Basically in the blue areas.

The square is located at (-1|-1),(1|-1),(1|1),(-1|1).
The circle has a radius of r = sqrt(2) and is centered at (0|0).

I already have scripts to:

• generate a random point on a circle (uniformly):

  float a = 2 * MathUtils.PI * MathUtils.random(1f); // angle between 0 and 2pi
  float r = radius * Math.sqrt(MathUtils.random(0, 1f)
  float x = r * MathUtils.cos(a);
  float y = r * MathUtils.sin(a);
  

• calculating the radius for a given angle to form a square:

  float r = (1/Math.sqrt(2)) / MathUtils.cos(((a+45)%90-45)/180*MathUtils.PI);
  

  with (1/Math.sqrt(2)) being half the side length of the square

Before anyone asks: I know that I could just re-generate points which are inside the green square until I get one that is outside, but I don't want to do it this way.

I appreciate any help. Thank you :)

Answer 1

It is rather hard to generate points only in region of sphere outside the cube (caps and wedges), so rejecting method looks reasonable.

But you can diminish the number of useless points, generating points in the ring only in 2D case and in spherical shell in 3D case.

So pseudocode might look as

 //2d
 SquaredR  = RandomUniformInRange(0.5, 1)
 R = Sqrt(SquaredR)

 //3d
 CubedR  = RandomUniformInRange(Pow(3, -3/2), 1)
 R = Pow(CubedR, 1/3)

 //generate point on the circle or on the sphere with radius R

 if Abs(x) > Sqrt(2)/2 or Sqrt(3)/3  and so on - reject

Having R, you can generate point on the sphere using any approach from here

Answer 2

Here is rough sketch of the idea. You select one quadrant to sample, say, one on the right.

First, sample angles from -pi/4 to pi/4

float a = -MathUtils.PI/4.0f + MathUtils.PI/2.0 * MathUtils.random(0.f,1.f);
float c = MathUtils.cos(a);
float s = MathUtils.sin(a);

Second, find minimum radius. With ray going from (0,0) at angle a will intersect quadrant line at x=1 at minimum

float rmin = 1.0f / c;
float rmax = Math.sqrt(2.0f);

Sample from rmin to rmax = sqrt(2), taking into account that for plane you sample squared radius and then use sqrt(), and for 3d space you sample cubed radius and then use cbrt().

float r2 = rmin*rmin + (rmax*rmax-rmin*rmin)*MathUtils.random(0.f,1.f);
float r  = Math.sqrt(r);

float x = r * c;
float y = r * s;

Now, we constructed (x,y) is a such way it is guaranteed to be in the right quadrant below circle and on the right of the x=1 line.

To cover all four quadrants just sample to which quadrant you will move the point

float q = MathUtils.random(0.f,1.f);
if (q < 0.25f) // top quadrant
    return (y, x);
if (q < 0.5f)  // left quadrant
    return (-x, y);
if (q < 0.75f)  // bottom quadrant
    return (y, -x);
return (x,y); // right quadrant

Please bear with me - my Java is quite rusty, and I have no ways to test the code.

In 3D case you'll have to deal with two angles, cubed radius, eight octants instead of four quadrants, but general logic is the same

UPDATE

I was wrong, sampling like I propose would lead to non-uniform point distribution.

From PDF:

PDF(phi, r) = S_(-\pi/4)^\phi d\phi S_1/cos(\phi)^\sqrt(2) r dr

One could get that we have to make \phi sampling non-uniform. Unfortunately, from U(0,1) to get to sampled \phi requires solving non-linear equation

\pi/2 (0.5*(\phi/\pi/4 + 1) - U(0,1)) = 0.5*(tan(phi) + 1) - U(0,1)

So algorithm would be:

• Sample U(0,1)
• Find appropriate \phi by solving equation above
• Find lower R boundary
• Sample R

Quick code (in Python, sorry) to plot this non-linear function

import numpy as np
import matplotlib.pyplot as plt

def f(phi, ksi):
    c   = 0.5 * np.pi
    c_2 = 0.5 * c

    left = c * (0.5 * (phi/c_2     + 1.0) - ksi)
    rght =     (0.5 * (np.tan(phi) + 1.0) - ksi)

    return left - rght

nof_points = 41
phi   = np.linspace(-0.25*np.pi, 0.25*np.pi, nof_points)
y0_00 = f(phi, 0.00)
y0_25 = f(phi, 0.25)
y0_50 = f(phi, 0.50)
y0_75 = f(phi, 0.75)
y0_99 = f(phi, 1.00)

plt.plot(phi, y0_00, 'ro', phi, y0_25, 'b+', phi, y0_50, 'gx', phi, y0_75, 'm.', phi, y0_99, 'y^')
plt.show()

and plotted functions for five values of U(0,1) (ksi in the code)

Sampling could be rearranged such that r sampling is non-linear, but it exhibits the same problem - need to solve non-linear equation with polynomial and trigonometric parts

UPDATE II

And just for the record, if you want to sample r first, then it has to be sampled from the solution of the non-linear equation: r² sec^-1(r) - sqrt(r² - 1) = U(0,1)*(\pi/2 - 1)

in the interval [1...sqrt(2)]

After solving it and finding sampled r, \phi could be sampled uniformly in the interval allowed by r: [-cos^-1(1/r) ... +cos^-1(1/r)]