Randomly sampling points within an octagon using Python

Question

How do I randomly sample 2d points uniformly from within an octagon using Python / Numpy? We can say that the octagon is centered at the origin (0, 0). The following is what I've done:

import numpy as np
import matplotlib.pyplot as plt

def sample_within_octagon(num_points):
    points = np.zeros((num_points, 2))

    # Generate random angle in radians
    angles = np.random.uniform(0, 2 * np.pi, size=(num_points,))

    # Calculate the maximum radius for the given angle
    # This is wrong.
    max_radii = 1.0 / np.sqrt(2) / np.cos(np.pi / 8 - angles % (np.pi / 4))

    # Generate random radius within the bounds of the octagon
    # Use square-root to prevent it from being more dense in center.
    radii = np.sqrt(np.random.uniform(0, max_radii))

    # Convert polar coordinates to Cartesian coordinates
    x = radii * np.cos(angles)
    y = radii * np.sin(angles)
    points[:, 0] = x
    points[:, 1] = y
    
    return points

num_points = 10000
random_points = sample_within_octagon(num_points)
plt.scatter(
    np.array(random_points)[:, 0], 
    np.array(random_points)[:, 1], s=1);
plt.axis('equal');

The above code is mostly correct, but the max_radii calculation is incorrect, because the edges are slightly curved outward.

I am not necessarily committed to the overall approach of the above algorithm, so any algorithm will do. Having said that, I would slightly prefer an approach that (like the above, if it had actually worked correctly) would generalize to 16-gons and so on.

Answer 1

In your code, the formula for max_radii needs a little modification, try the following:

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate


def sample_within_octagon(num_points, inv_transform_evals=10000):
  points = np.zeros((num_points, 2))
  # Angle offset for each octagon segment
  offset = np.pi / 8.0
  # Generate random angle in radians
  max_radii_in = np.linspace(0, 2 * np.pi, inv_transform_evals)
  max_radii_out = 1 / np.cos(np.abs(max_radii_in % (np.pi / 4) - offset))
  max_radii_cdf = np.cumsum(max_radii_out / max_radii_out.sum())
  f = interpolate.interp1d(np.array([0.] + list(max_radii_cdf)),
                           np.array([0.] + list(max_radii_in)))
  angles_out = np.random.uniform(0, 1, num_points)
  angles = f(angles_out)
  # Calculate max radius based on octagon geometry
  max_radii = 1 / np.cos(np.abs(angles % (np.pi / 4) - offset))
  # Generate random radius with square root scaling
  radii = np.sqrt(np.random.uniform(0, 1, num_points)) * max_radii
  # Convert to Cartesian coordinates
  points[:, 0] = radii * np.cos(angles)
  points[:, 1] = radii * np.sin(angles)
  return points


# Generate and plot points
num_points = 10000
points = sample_within_octagon(num_points)
plt.scatter(points[:, 0], points[:, 1], s=1)
plt.axis('equal')
plt.show()

Note: The above solution has been modified by the OP - @calmcc based on suggestions in the comments of the question.

Answer 2

I would like to propose alternative, which is based on triangulation of the octagon. You split it into 8 triangles, choice triangle with equal probabilities, and then sample from known method of uniform sampling in triangle, see Generate random locations within a triangular domain for details.

Could be used then for hexagon, 16gon and in general for any triangulated area (well, triangle choice would be with probability proportional to area, but that would be the only modification)

Code, Windows x64 Python 3.11

import numpy as np
import matplotlib.pyplot as plt
    
X = 0
Y = 1

rng = np.random.default_rng(135797531)

def trisample(A, B, C): # https://stackoverflow.com/questions/47410054/generate-random-locations-within-a-triangular-domain/47425047#47425047

    """
    Given three vertices A, B, C, 
    sample point uniformly in the triangle
    """
    r1 = rng.random()
    r2 = rng.random()

    s1 = np.sqrt(r1)

    x = A[X] * (1.0 - s1) + B[X] * (1.0 - r2) * s1 + C[X] * r2 * s1
    y = A[Y] * (1.0 - s1) + B[Y] * (1.0 - r2) * s1 + C[Y] * r2 * s1

    return (x, y)

R = 10
M = 8 # need octagone

A = np.empty((M,3))
B = np.empty((M,3))
C = np.empty((M,3))

delta_phi = 2.0 * np.pi / M

for k in range(0, M): # could be vectorized but we do it only once
    phil = k * delta_phi
    phir = (k+1) * delta_phi

    A[k, X] = 0.0
    A[k, Y] = 0.0
    
    B[k, X] = R*np.cos(phil)
    B[k, Y] = R*np.sin(phil)
    
    C[k, X] = R*np.cos(phir)
    C[k, Y] = R*np.sin(phir)
    

def sample_within_octagon(num_points):
    trindices = rng.choice(M, size=num_points, replace=True)
    
    points = np.empty((num_points, 2))
    
    for k in range(0, num_points): # could be vectorized
        idx = trindices[k]
        points[k,X], points[k,Y] = trisample(A[idx,:], B[idx,:], C[idx,:])

    return points

num_points = 10000
random_points = sample_within_octagon(num_points)
plt.scatter(
    np.array(random_points)[:, 0], 
    np.array(random_points)[:, 1], s=1);
plt.axis('equal');

and it will produce nice picture like below