使用Python在八边形内随机采样点

Question

如何使用 Python / Numpy 从八边形内均匀随机采样 2d 点？我们可以说八边形的中心位于原点 (0, 0)。以下是我所做的：

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');

上面的代码大部分是正确的，但是 max_radii 计算不正确，因为边缘稍微向外弯曲。

我不一定致力于上述算法的整体方法，所以任何算法都可以。话虽如此，我稍微更喜欢一种可以推广到 16 边形等的方法（如上面的方法，如果它实际上工作正常的话）。

Answer 1

在您的代码中，

max_radii

的公式不正确，并且没有必要使用

np.sqrt

作为半径。尝试以下操作：

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 normalized angle within each sector
  norm_angles = np.abs((angles % (np.pi / 4)) - np.pi / 8)
  # Calculate the maximum radius for the given angle
  max_radii = 1 / (np.cos(norm_angles) + np.sin(norm_angles))
  # Generate random radius within the bounds of the octagon
  radii = np.random.uniform(0, max_radii)
  # Convert polar coordinates to Cartesian coordinates
  points[:, 0] = radii * np.cos(angles)
  points[:, 1] = radii * np.sin(angles)
  return points


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