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

Question

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

3

如何使用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边形等等。

- calmcc

1

有两个问题。第一个是半径计算 - 首先尝试编写和调试一个单独处理单个输入角度的函数，以确保您具有正确的基本逻辑。第二个问题是，“选择一个随机角度，然后从中心选择一个随机距离”不会得到均匀的结果，即使进行了平方根修正 - 因为“半径”（实际上不是一个圆的半径）根据角度而变化，所以随机角度选择会对沿着较短“半径”的点产生偏向。 - Karl Knechtel

1

一个常见的做法是首先创建一个包围目标形状的矩形，然后在矩形内随机生成点，最后筛选出在目标形状外部的点。假设这不是一个数学练习，你需要多精确呢？ - ken

@ken 那种方法（称为拒绝抽样）是完全准确的；这不是问题所在。问题是它不是一个确定性算法：假设输入是真正随机的，它不能保证在有限时间内终止。 - Karl Knechtel

2

@KarlKnechtel 我已经澄清了问题，表示我会接受任何有效的答案，尽管我对我尝试的方法有一点情感依恋。 - calmcc

1

我为你在Codidact上问了这个数学问题：https://math.codidact.com/posts/289536 - Karl Knechtel

显示剩余5条评论

3个回答

2

我想提出一种基于八边形三角化的替代方案。您将其分割成8个三角形，以等概率选择一个三角形，然后从已知的均匀采样方法中进行采样，详细信息请参见在三角形域内生成随机位置。

然后可以用于六边形、16边形，以及一般的任何三角化区域（嗯，三角形的选择将与面积成比例的概率有关，但这将是唯一的修改）。

代码，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://dev59.com/5qbja4cB1Zd3GeqPn-76#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');

并且它会生成像下面这样的漂亮图片

- Severin Pappadeux

这在我看来简单得多，并且符合我在 Codidact 上尝试询问时得出的相同结论。 - Karl Knechtel

@KarlKnechtel 它可能会更快，尤其是如果您不对三角形索引进行采样，而只对每个三角形中的总点数的1/8进行采样。采样例程只有一个sqrt()，所有sin()/cos()调用都在准备阶段完成。 - Severin Pappadeux

0

为了在八边形（或者任何多边形）内均匀采样二维点，我想你必须从一开始就使用X和Y坐标。

-计算多边形的边界框 -在边界框内生成均匀分布的随机点 -仅保留多边形内的点 -确保返回正确数量的点（迭代直到达到所需点的数量）

- rellachs

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Sash Sinha · Accepted Answer

在你的代码中，max_radii 的公式需要稍作修改，请尝试以下方法：

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()

注意：上述解决方案已由提问者@calmcc根据问题评论中的建议进行了修改。