Generate a random point within a circle (uniformly

Question

I need to generate a uniformly random point within a circle of radius R.

I realize that by just picking a uniformly random angle in the interval [0 ... 2π), and uniformly random radius in the interval (0 ... R) I would end up with more points towards the center, since for two given radii, the points in the smaller radius will be closer to each other than for the points in the larger radius.

I found a blog entry on this over here but I don't understand his reasoning. I suppose it is correct, but I would really like to understand from where he gets (2/R²)×r and how he derives the final solution.

---

Update: 7 years after posting this question I still hadn't received a satisfactory answer on the actual question regarding the math behind the square root algorithm. So I spent a day writing an answer myself. Link to my answer.

Answer 1

The reason why the naive solution doesn't work is that it gives a higher probability density to the points closer to the circle center. In other words the circle that has radius r/2 has probability r/2 of getting a point selected in it, but it has area (number of points) pi*r^2/4.

Therefore we want a radius probability density to have the following property:

The probability of choosing a radius smaller or equal to a given r has to be proportional to the area of the circle with radius r. (because we want to have a uniform distribution on the points and larger areas mean more points)

In other words we want the probability of choosing a radius between [0,r] to be equal to its share of the overall area of the circle. The total circle area is pi*R^2, and the area of the circle with radius r is pi*r^2. Thus we would like the probability of choosing a radius between [0,r] to be (pi*r^2)/(pi*R^2) = r^2/R^2.

Now comes the math:

The probability of choosing a radius between [0,r] is the integral of p(r) dr from 0 to r (that's just because we add all the probabilities of the smaller radii). Thus we want integral(p(r)dr) = r^2/R^2. We can clearly see that R^2 is a constant, so all we need to do is figure out which p(r), when integrated would give us something like r^2. The answer is clearly r * constant. integral(r * constant dr) = r^2/2 * constant. This has to be equal to r^2/R^2, therefore constant = 2/R^2. Thus you have the probability distribution p(r) = r * 2/R^2

Note: Another more intuitive way to think about the problem is to imagine that you are trying to give each circle of radius r a probability density equal to the proportion of the number of points it has on its circumference. Thus a circle which has radius r will have 2 * pi * r "points" on its circumference. The total number of points is pi * R^2. Thus you should give the circle r a probability equal to (2 * pi * r) / (pi * R^2) = 2 * r/R^2. This is much easier to understand and more intuitive, but it's not quite as mathematically sound.

Answer 2

I used once this method: This may be totally unoptimized (ie it uses an array of point so its unusable for big circles) but gives random distribution enough. You could skip the creation of the matrix and draw directly if you wish to. The method is to randomize all points in a rectangle that fall inside the circle.

bool[,] getMatrix(System.Drawing.Rectangle r) {
    bool[,] matrix = new bool[r.Width, r.Height];
    return matrix;
}

void fillMatrix(ref bool[,] matrix, Vector center) {
    double radius = center.X;
    Random r = new Random();
    for (int y = 0; y < matrix.GetLength(0); y++) {
        for (int x = 0; x < matrix.GetLength(1); x++)
        {
            double distance = (center - new Vector(x, y)).Length;
            if (distance < radius) {
                matrix[x, y] = r.NextDouble() > 0.5;
            }
        }
    }

}

private void drawMatrix(Vector centerPoint, double radius, bool[,] matrix) {
    var g = this.CreateGraphics();

    Bitmap pixel = new Bitmap(1,1);
    pixel.SetPixel(0, 0, Color.Black);

    for (int y = 0; y < matrix.GetLength(0); y++)
    {
        for (int x = 0; x < matrix.GetLength(1); x++)
        {
            if (matrix[x, y]) {
                g.DrawImage(pixel, new PointF((float)(centerPoint.X - radius + x), (float)(centerPoint.Y - radius + y)));
            }
        }
    }

    g.Dispose();
}

private void button1_Click(object sender, EventArgs e)
{
    System.Drawing.Rectangle r = new System.Drawing.Rectangle(100,100,200,200);
    double radius = r.Width / 2;
    Vector center = new Vector(r.Left + radius, r.Top + radius);
    Vector normalizedCenter = new Vector(radius, radius);
    bool[,] matrix = getMatrix(r);
    fillMatrix(ref matrix, normalizedCenter);
    drawMatrix(center, radius, matrix);
}

Answer 3

Let ρ (radius) and φ (azimuth) be two random variables corresponding to polar coordinates of an arbitrary point inside the circle. If the points are uniformly distributed then what is the disribution function of ρ and φ?

For any r: 0 < r < R the probability of radius coordinate ρ to be less then r is

P[ρ < r] = P[point is within a circle of radius r] = S1 / S0 =(r/R)²

Where S1 and S0 are the areas of circle of radius r and R respectively. So the CDF can be given as:

          0          if r<=0
  CDF =   (r/R)**2   if 0 < r <= R
          1          if r > R

And PDF:

PDF = d/dr(CDF) = 2 * (r/R**2) (0 < r <= R).

Note that for R=1 random variable sqrt(X) where X is uniform on [0, 1) has this exact CDF (because P[sqrt(X) < y] = P[x < y**2] = y**2 for 0 < y <= 1).

The distribution of φ is obviously uniform from 0 to 2*π. Now you can create random polar coordinates and convert them to Cartesian using trigonometric equations:

x = ρ * cos(φ)
y = ρ * sin(φ)

Can't resist to post python code for R=1.

from matplotlib import pyplot as plt
import numpy as np

rho = np.sqrt(np.random.uniform(0, 1, 5000))
phi = np.random.uniform(0, 2*np.pi, 5000)

x = rho * np.cos(phi)
y = rho * np.sin(phi)

plt.scatter(x, y, s = 4)

You will get

Answer 4

I am still not sure about the exact '(2/R2)×r' but what is apparent is the number of points required to be distributed in given unit 'dr' i.e. increase in r will be proportional to r2 and not r.

check this way...number of points at some angle theta and between r (0.1r to 0.2r) i.e. fraction of the r and number of points between r (0.6r to 0.7r) would be equal if you use standard generation, since the difference is only 0.1r between two intervals. but since area covered between points (0.6r to 0.7r) will be much larger than area covered between 0.1r to 0.2r, the equal number of points will be sparsely spaced in larger area, this I assume you already know, So the function to generate the random points must not be linear but quadratic, (since number of points required to be distributed in given unit 'dr' i.e. increase in r will be proportional to r2 and not r), so in this case it will be inverse of quadratic, since the delta we have (0.1r) in both intervals must be square of some function so it can act as seed value for linear generation of points (since afterwords, this seed is used linearly in sin and cos function), so we know, dr must be quadratic value and to make this seed quadratic, we need to originate this values from square root of r not r itself, I hope this makes it little more clear.

Answer 5

How to generate a random point within a circle of radius R:

r = R * sqrt(random())
theta = random() * 2 * PI

(Assuming random() gives a value between 0 and 1 uniformly)

If you want to convert this to Cartesian coordinates, you can do

x = centerX + r * cos(theta)
y = centerY + r * sin(theta)

Why sqrt(random())?

Let's look at the math that leads up to sqrt(random()). Assume for simplicity that we're working with the unit circle, i.e. R = 1.

The average distance between points should be the same regardless of how far from the center we look. This means for example, that looking on the perimeter of a circle with circumference 2 we should find twice as many points as the number of points on the perimeter of a circle with circumference 1.

                

Since the circumference of a circle (2πr) grows linearly with r, it follows that the number of random points should grow linearly with r. In other words, the desired probability density function (PDF) grows linearly. Since a PDF should have an area equal to 1 and the maximum radius is 1, we have

                

So we know how the desired density of our random values should look like. Now: How do we generate such a random value when all we have is a uniform random value between 0 and 1?

We use a trick called inverse transform sampling

1. From the PDF, create the cumulative distribution function (CDF)
2. Mirror this along y = x
3. Apply the resulting function to a uniform value between 0 and 1.

Sounds complicated? Let me insert a yellow box with a little sidetrack that conveys the intuition:

Suppose we want to generate a random point with the following distribution:

                

That is

• 1/5 of the points uniformly between 1 and 2, and
• 4/5 of the points uniformly between 2 and 3.

The CDF is, as the name suggests, the cumulative version of the PDF. Intuitively: While PDF(x) describes the number of random values at x, CDF(x) describes the number of random values less than x.

In this case the CDF would look like:

                

To see how this is useful, imagine that we shoot bullets from left to right at uniformly distributed heights. As the bullets hit the line, they drop down to the ground:

                

See how the density of the bullets on the ground correspond to our desired distribution! We're almost there!

The problem is that for this function, the y axis is the output and the x axis is the input. We can only "shoot bullets from the ground straight up"! We need the inverse function!

This is why we mirror the whole thing; x becomes y and y becomes x:

                

We call this CDF^-1. To get values according to the desired distribution, we use CDF^-1(random()).

…so, back to generating random radius values where our PDF equals 2x.

Step 1: Create the CDF:

Since we're working with reals, the CDF is expressed as the integral of the PDF.

CDF(x) = ∫ 2x = x²

Step 2: Mirror the CDF along y = x:

Mathematically this boils down to swapping x and y and solving for y:

CDF:     y = x²
Swap:   x = y²
Solve:   y = √x
CDF^-1:  y = √x

Step 3: Apply the resulting function to a uniform value between 0 and 1

CDF^-1(random()) = √random()

Which is what we set out to derive :-)

Answer 6

First we generate a cdf[x] which is

The probability that a point is less than distance x from the centre of the circle. Assume the circle has a radius of R.

obviously if x is zero then cdf[0] = 0

obviously if x is R then the cdf[R] = 1

obviously if x = r then the cdf[r] = (Pi r^2)/(Pi R^2)

This is because each "small area" on the circle has the same probability of being picked, So the probability is proportionally to the area in question. And the area given a distance x from the centre of the circle is Pi r^2

so cdf[x] = x^2/R^2 because the Pi cancel each other out

we have cdf[x]=x^2/R^2 where x goes from 0 to R

So we solve for x

R^2 cdf[x] = x^2

x = R Sqrt[ cdf[x] ]

We can now replace cdf with a random number from 0 to 1

x = R Sqrt[ RandomReal[{0,1}] ]

Finally

r = R Sqrt[  RandomReal[{0,1}] ];
theta = 360 deg * RandomReal[{0,1}];
{r,theta}

we get the polar coordinates {0.601168 R, 311.915 deg}