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.

You've probably thought about solutions to this problem. What are they, how do you think that they're failing, what code did you try to write and how exactly does it fail to do what you expect? Without that, it looks a lot like a "do my job for me for free" question... — Thierry Lathuille, Aug 26 '23 at 13:50
@KarlKnechtel -- I've showed where I'm stuck now. Fwiw, I'm not a student -- this is for a personal project. — calmcc, Aug 26 '23 at 14:15
There are two problems. The first is the radius calculation - try first writing and debugging a separate function that works with a single input angle, to make sure you have correct underlying logic. The second is that "choose a random angle, then a random distance from the center" won't give a uniform result, even after the square-root correction - because the "radius" (not really a radius because it's not a circle) varies according to the angle, so the random-angle selection creates a bias towards points that along the shorter "radii". — Karl Knechtel, Aug 26 '23 at 14:25
A common practice is to first create a rectangle that encloses the target shape, then randomly generate points within the rectangle, and finally filter out points that are outside the target shape. Assuming this is not a math exercise, how accurate do you need to be? — ken, Aug 26 '23 at 14:27
(If the question is specifically about the *algorithm*, then I don't think we actually need code. But you should be clear about whether you want to debug this approach, or whether you want a different approach laid out.) — Karl Knechtel, Aug 26 '23 at 14:27
@ken that approach (called [*rejection sampling*](https://en.wikipedia.org/wiki/Rejection_sampling)) is perfectly accurate; that isn't the concern. The concern is that it is not a deterministic algorithm: assuming a true random input, it is not guaranteed to terminate in finite time. — Karl Knechtel, Aug 26 '23 at 14:29
@KarlKnechtel I've clarified the question to say that I'll accept any answer that works, even though I do have a bit of an attachment to the approach that I was trying. — calmcc, Aug 26 '23 at 14:38
... Although I think the algorithm question is better suited to [math.se]. — Karl Knechtel, Aug 26 '23 at 14:39
I asked the Math question for you on Codidact: https://math.codidact.com/posts/289536 — Karl Knechtel, Aug 26 '23 at 20:39

Sash Sinha · Accepted Answer · 2023-08-26T16:09:06.220

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.

Really cool, thanks! I've edited your solution, using inverse transform sampling, to address KarlKnechtel's observation that random-angle selection creates a bias towards points that along the shorter "radii". — calmcc, Aug 26 '23 at 15:36

score 2 · Answer 2 · answered Aug 26 '23 at 15:33

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

This is much simpler IMO, and follows from the same conclusion that was reached when I tried asking on Codidact. — Karl Knechtel, Aug 26 '23 at 20:41
@KarlKnechtel It would be probably faster as well especially if you don't sample triangle index and just sample 1/8 of the total points in each triangle. Sampling routine is one sqrt(), all sin()/cos() calls are done at preparation stage — Severin Pappadeux, Aug 26 '23 at 21:07

score 0 · Answer 3 · answered Aug 26 '23 at 14:55

In order to uniformly sample 2D points within an octagon (or any polygon for that matter), I guess you have to stick with X and Y coordinates right from the start.

calculate the bounding box of the polygon
generate uniformly distributed random points within the bounding box
keep only the points within the polygon
make sure you return the correct number of points (iterate until you reach the desired number of points)

Randomly sampling points within an octagon using Python

3 Answers3