import matplotlib.pyplot as plt
import numpy as np
# Number of points to sample
= 1000
num_points
# Arrays to store the sampled points
= []
x_points = []
y_points
for _ in range(num_points):
= np.random.uniform(0, 1)
r = np.random.uniform(0, r)
x = np.sqrt(r**2 - x**2)
y
x_points.append(x)
y_points.append(y)
# Plotting the points
=(6, 6))
plt.figure(figsize='turquoise', s=2)
plt.scatter(x_points, y_points, color0, 1)
plt.xlim(0, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Initial Sampling - Points Clustered Around Center')
plt.title(True)
plt.grid( plt.show()
1 Introduction
2 Introduction
YouTube recommendations can sometimes lead to surprisingly insightful videos. One such example is the video The BEST Way to Find a Random Point in a Circle by nubDotDev, which was originally submitted to 3blue1brown’s Summer of Math Exposition. This video addresses an interesting problem: How would you write an algorithm to sample points uniformly within a circle?
Before reading this post, I recommend watching the video, as it covers rejection sampling and polar coordinate methods to solve this problem. In the following sections, I will present another solution that neither relies on rejection sampling nor polar coordinates, but instead utilizes the equation of a circle directly. Why would you want to know about all this? I’m not sure, but maybe this simple problem stimulates your curiosity like it did for me, so I recommend doing it for this purpose alone.
As already mentioned, there exist many solutions to this problem, like the rejection-sampling method. This post is about the result of my own ideas of thinking about the problem, but of course, these are not new and can be found in the existing approaches.
3 Uniform Sampling
3.1 Naive Approach
The equation for a circle of radius \(r\) centered at the origin is given by:
\[ y^2 + x^2 = r^2 \tag{1}\]
Using Equation 1, random points on a circle can be sampled by selecting random values for \(x \in [0,r]\) and calculating the corresponding \(y\) value (or vice versa). This only gives us points lying on the circle, whereas our objective is to sample uniformly from the entire area within the circle. This can be achieved by additionally sampling the radius as \(r_{\text{sampled}} \in [0,r]\). To summarize, the following procedure can be followed as our first attempt to sample points uniformly within a circle of radius \(r = 1\):
- Sample \(r_{\text{sampled}} \in [0,1]\)
- Sample \(x \in [0, r_{\text{sampled}}]\)
- Calculate \(y\) (Equation 1)
- Plot the point \((x,y)\)
To keep the code simple, we’ll just create a quarter circle, that is use \(x,y \in [0,1]\).
3.2 Circumference sampling
This method is simple, but it has a couple of problems, resulting in a non-uniform density of points. The most obvious problem is, that the points are denser around the circles’ center. The problem arises because the circumference of the circle grows with \(r_{\text{sampled}}\) but the number of samples \(k\) doesn’t, as shown in this image.
To fix this, we can simply generate more points as the circle grows. Specifically, we’ll use the circumference \(c\) as a scaling factor.
\[ c = 2 \pi r_{\text{sampled}} \tag{2}\]
def points_based_on_circumference(r, scale=1):
return int(np.ceil(2 * np.pi * r * scale))
= []
x_points = []
y_points = 1000
num_radii
for _ in range(num_radii):
= np.random.uniform(0, 1)
r = points_based_on_circumference(r)
num_points for _ in range(num_points):
= np.random.uniform(0, r)
x = np.sqrt(r**2 - x**2)
y
x_points.append(x)
y_points.append(y)
=(6, 6))
plt.figure(figsize='turquoise', s=2)
plt.scatter(x_points, y_points, color0, 1)
plt.xlim(0, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Improved Distribution with Circumference-Based Sampling')
plt.title(True)
plt.grid( plt.show()
3.3 Inverse
Now we notice the second problem: It seems like there are fewer points at the bottom of the diagram. To explain what’s going on, it’s easier to look at a single circle generated with this approach.
Code
# Number of points to sample
= 100
num_points
# Arrays to store the sampled points
= []
x_points = []
y_points
for _ in range(num_points):
# Sample r from [0, 1]
= 1
r
# Sample x from [0, r]
= np.random.uniform(0, r)
x
# Calculate y based on the circle formula X^2 + Y^2 = r^2
= np.sqrt(r**2 - x**2)
y
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='turquoise', s=5)
plt.scatter(x_points, y_points, color
# Set the limits and labels
0, 1)
plt.xlim(0, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Points on a Circle')
plt.title(True)
plt.grid(
# Show the plot
plt.show()
Because the circle is drawn by calculating \(y\) via \(x\), which is uniformly sampled, there are simply more values on the “roof” of the circle than on the side, because from the perspective of \(x\), the area is wider. We can fix t his problem by also considering the perspective of \(y\), that is inverting \(x\) and \(y\) by chance.
import random
# Number of points to sample
= 100
num_points
# Arrays to store the sampled points
= []
x_points = []
y_points
for _ in range(num_points):
# Sample r from [0, 1]
= 1
r
if random.random() < 0.5:
# Sample x from [0, r]
= np.random.uniform(0, r)
x
# Calculate y based on the circle formula X^2 + Y^2 = r^2
= np.sqrt(r**2 - x**2)
y else:
= np.random.uniform(0, r)
y = np.sqrt(r**2 - y**2)
x
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='turquoise', s=5)
plt.scatter(x_points, y_points, color
# Set the limits and labels
0, 1)
plt.xlim(0, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Points on a Circle')
plt.title(True)
plt.grid(
# Show the plot
plt.show()
3.4 Full circle
Using the \(c\) as a factor for the number of samples \(k\) and by inversing \(x\) and \(y\) by chance, we can finally generate a full circle. This is done by sampling \(x \in [-r_{\text{sampled}}, r_{\text{sampled}}]\) and mirroring the point on the x-axis by chance.
import matplotlib.pyplot as plt
import numpy as np
# Function to determine the number of points based on the circumference
def points_based_on_circumference(r, scale=1):
return int(np.ceil(2 * np.pi * r * scale))
# Arrays to store the sampled points
= []
x_points = []
y_points
# Number of different radii to sample
= 1000
num_radii
for _ in range(num_radii):
# Sample r from [0, 1]
= np.random.uniform(0, 1)
r
# Number of points to sample for this radius
= points_based_on_circumference(r)
num_points
for _ in range(num_points):
if random.random() < 0.5:
# Sample x from [0, r]
= np.random.uniform(-r, r)
x
# Calculate y based on the circle formula X^2 + Y^2 = r^2
= np.sqrt(r**2 - x**2)
y if np.random.rand() > 0.5:
= -y
y else:
= np.random.uniform(-r, r)
y = np.sqrt(r**2 - y**2)
x
if np.random.rand() > 0.5:
= -x
x
# Randomly decide the sign of y to cover both the top and bottom halves
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='turquoise', s=2)
plt.scatter(x_points, y_points, color
# Set the limits and labels
-1, 1)
plt.xlim(-1, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Points on Circles with Different Radii')
plt.title(True)
plt.grid(
# Show the plot
plt.show()
With a bit more refactoring, we can also make the code more compact and reusable.
Refactored code
import matplotlib.pyplot as plt
import numpy as np
import math
= 10
target_r
= 500
r_granularity = 1
c_granularity
# Function to determine the number of points based on the circumference
def points_based_on_circumference(r, granularity):
return int(np.ceil(2 * np.pi * r * granularity))
# Arrays to store the sampled points
= []
x_points = []
y_points
# Number of different radii to sample
for _ in range(r_granularity):
= np.random.uniform(0, target_r)
r
# Number of points to sample for this radius
= points_based_on_circumference(r, c_granularity)
num_points
for _ in range(num_points):
# Sample and calculate vars
= np.random.uniform(-r, r)
dep = np.sqrt(r**2 - dep**2)
indep
# Mirror on x-axis
if np.random.rand() > 0.5:
= -indep
indep
# Assign indep. and dep. var
= (dep, indep) if np.random.rand() > 0.5 else (indep, dep)
x, y
# Append the point to the arrays
x_points.append(x) y_points.append(y)
3.5 Blurring the lines
However, with this approach, we will encounter a final problem. For small values of r_granularity
, it becomes clear that the points are sampled on separate circles.
Code
= 10
target_r
= 25
r_granularity = 5
c_granularity
# Function to determine the number of points based on the circumference
def points_based_on_circumference(r, granularity):
return int(np.ceil(2 * np.pi * r * granularity))
# Arrays to store the sampled points
= []
x_points = []
y_points
# Number of different radii to sample
for _ in range(r_granularity):
= np.random.uniform(0, target_r)
r
# Number of points to sample for this radius
= points_based_on_circumference(r, c_granularity)
num_points
for _ in range(num_points):
# Sample and calculate vars
= np.random.uniform(-r, r)
dep = np.sqrt(r**2 - dep**2)
indep
# Mirror on x-axis
if np.random.rand() > 0.5:
= -indep
indep
# Assign indep. and dep. var
= (dep, indep) if np.random.rand() > 0.5 else (indep, dep)
x, y
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='turquoise', s=2)
plt.scatter(x_points, y_points, color
# Set the limits and labels
-target_r, target_r)
plt.xlim(-target_r, target_r)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('target_r = 10, r_granularity = 25, c_granularity = 5')
plt.title(True)
plt.grid(
# Show the plot
plt.show()
To solve this problem, we can replace the scaling factor, which is calculated using the circumference \(c\) with a probability. So instead of deterministically drawing \(k = c_{\text{sampled}} = 2 \pi r_{\text{sampled}}\) points on a circle with radius \(r_{\text{sampled}}\) we will only draw one point with probability \(p = \frac{c_{\text{sampled}}}{c}\).
= 1
target_r
= target_r * 10000
r_granularity
# Function to determine the number of points based on the circumference
def p_acc_based_on_circumference(r):
return (2 * np.pi * r) / (2 * np.pi * target_r)
# Arrays to store the sampled points
= []
x_points = []
y_points
# Number of different radii to sample
for _ in range(r_granularity):
= np.random.uniform(0, target_r)
r
# Number of points to sample for this radius
= p_acc_based_on_circumference(r)
p_acc
if np.random.rand() < p_acc:
# Sample and calculate vars
= np.random.uniform(-r, r)
dep = np.sqrt(r**2 - dep**2)
indep
# Mirror on x-axis
if np.random.rand() > 0.5:
= -indep
indep
# Assign indep. and dep. var
= (dep, indep) if np.random.rand() > 0.5 else (indep, dep)
x, y
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='turquoise', s=2)
plt.scatter(x_points, y_points, color
# Set the limits and labels
-target_r, target_r)
plt.xlim(-target_r, target_r)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Final algorithm')
plt.title(True)
plt.grid(
# Show the plot
plt.show()
4 Extra
The approach introduced in the previous sections works not only for circles but also for other shapes. The different components like Equation 1 or the ways of mirroring can be adjusted to create entirely new shapes with uniformly sampled points. I will leave it up to the reader to find the changes made to each of these following shapes, resulting in this variety.
Code
import matplotlib.pyplot as plt
import numpy as np
# Function to determine the number of points based on the circumference
def points_based_on_circumference(r, scale=100):
return int(np.ceil(2 * np.pi * r * scale))
# Arrays to store the sampled points
= []
x_points = []
y_points
# Number of different radii to sample
= 1000
num_radii
for _ in range(num_radii):
# Sample r from [0, 1]
= np.random.uniform(0, 1)
r
# Number of points to sample for this radius
= points_based_on_circumference(r) // 20
num_points
for _ in range(num_points):
# Sample x from [-r, r]
= np.random.uniform(-r, r)
x
# Calculate y based on the circle formula X^2 + Y^2 = r^2
= np.sqrt(r - x**2)
y
# Randomly decide the sign of y to cover both the top and bottom halves
if np.random.rand() > 0.5:
= -y
y
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='firebrick', s=2)
plt.scatter(x_points, y_points, color
# Set the limits and labels
-1, 1)
plt.xlim(-1, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Modified shape 1')
plt.title(True)
plt.grid(
# Show the plot
plt.show()
Code
import matplotlib.pyplot as plt
import numpy as np
# Function to determine the number of points based on the circumference
def points_based_on_circumference(r, scale=100):
return int(np.ceil(2 * np.pi * r * scale))
# Arrays to store the sampled points
= []
x_points = []
y_points
# Number of different radii to sample
= 1000
num_radii
for _ in range(num_radii):
# Sample r from [0, 1]
= np.random.uniform(0, 1)
r
# Number of points to sample for this radius
= points_based_on_circumference(r) // 20
num_points
for _ in range(num_points):
# Sample x from [-r, r]
= np.random.uniform(-r, r)
x
# Calculate y based on the modified equation X + Y = r
= r - abs(x)
y
# Randomly decide the sign of y to cover both the top and bottom halves
if np.random.rand() > 0.5:
= -y
y
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='firebrick', s=2)
plt.scatter(x_points, y_points, color
# Set the limits and labels
-1, 1)
plt.xlim(-1, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Modified shape 2')
plt.title(True)
plt.grid(
# Show the plot
plt.show()
Code
import matplotlib.pyplot as plt
import numpy as np
# Function to determine the number of points based on the circumference
def points_based_on_circumference(r, scale=100):
return int(np.ceil(2 * np.pi * r * scale))
# Arrays to store the sampled points
= []
x_points = []
y_points
# Number of different radii to sample
= 5000
num_radii
for _ in range(num_radii):
# Sample r from [0, 1]
= np.random.uniform(0, 1)
r
# Number of points to sample for this radius
= points_based_on_circumference(r) // 100
num_points
for _ in range(num_points):
# Sample x from [-r, r]
= np.random.uniform(-r, r)
x
# Calculate y based on the equation x + y = r
= r - abs(x)
y
# Randomly decide the sign of y to cover both the top and bottom halves
if np.random.rand() > 0.5:
= -y
y
# Append the point to the arrays
x_points.append(x)
y_points.append(y)
# Create the plot
=(6, 6))
plt.figure(figsize='firebrick', s=2)
plt.scatter(x_points, y_points, color
# Set the limits and labels
-1, 1)
plt.xlim(-1, 1)
plt.ylim('X')
plt.xlabel('Y')
plt.ylabel('Modified shape 3')
plt.title(True)
plt.grid(
# Show the plot
plt.show()