Search code examples
pythonrandomshapes

Randomly sampling points within an octagon using Python

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.

octagon

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.

Solution

  • 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.

    octagon_v3