算法：计算椭圆内的伪随机点

Question

对于我正在制作的简单粒子系统，我需要给定一个具有宽度和高度的椭圆，计算位于该椭圆内的随机点 X, Y。

现在我的数学不是最好的，所以我想在这里问是否有人可以给我指出正确的方向。

也许正确的方法是在宽度范围内选择一个随机浮点，将其作为 X 并从中计算 Y 值？

Answer 1

1. 在半径为 1 的圆内生成一个随机点。这可以通过在区间

   phi

   中取随机角度

   [0, 2*pi)

   和在区间

   rho

   中随机值

   [0, 1)

   并计算
   来完成

   x = sqrt(rho) * cos(phi)
   y = sqrt(rho) * sin(phi)
   

   公式中的平方根确保圆内均匀分布。

2. 将

   x

   和

   y

   缩放到椭圆的尺寸

   x = x * width/2.0
   y = y * height/2.0
   

Answer 2

使用拒绝采样：在椭圆周围的矩形中选择一个随机点。通过检查 (x-x0)^2/a^2+(y-y0)^2/b^2-1 的符号来测试该点是否在椭圆内部。如果该点不在内部，则重复此操作。 （这假设椭圆与坐标轴对齐。类似的解决方案适用于一般情况，但当然更复杂。）

Answer 3

通过仔细考虑极坐标形式的定义，也可以在不使用拒绝采样的情况下生成椭圆内的点。从 wikipedia 椭圆的极坐标形式由

给出

直观地说，半径越大的地方，我们应该更频繁地采样极角θ。用更数学的方式来说，随机变量 θ 的 PDF 应该是 p(θ) dθ = dA / A，其中 dA 是角度 θ 且宽度为 dθ 的单个线段的面积。使用极角面积 dA = 1/2 r² dθ 的方程，椭圆面积为 π a b，则 PDF 变为

要从此 PDF 中随机采样，一种直接方法是逆 CDF 技术。这需要计算累积密度函数（CDF），然后对该函数求逆。使用 Wolfram Alpha 获得不定积分，然后将其求逆得到

的逆 CDF

其中 u 在 0 和 1 之间运行。因此，要对随机角度 θ 进行采样，您只需生成一个在 0 和 1 之间的均匀随机数 u，并将其代入该逆 CDF 方程。

要获得随机半径，可以使用适用于圆的相同技术（例如，请参见（均匀）在圆内生成随机点）。

这里是一些实现该算法的示例 Python 代码：

import numpy
import matplotlib.pyplot as plt
import random

# Returns theta in [-pi/2, 3pi/2]
def generate_theta(a, b):
    u = random.random() / 4.0
    theta = numpy.arctan(b/a * numpy.tan(2*numpy.pi*u))

    v = random.random()
    if v < 0.25:
        return theta
    elif v < 0.5:
        return numpy.pi - theta
    elif v < 0.75:
        return numpy.pi + theta
    else:
        return -theta

def radius(a, b, theta):
    return a * b / numpy.sqrt((b*numpy.cos(theta))**2 + (a*numpy.sin(theta))**2)

def random_point(a, b):
    random_theta = generate_theta(a, b)
    max_radius = radius(a, b, random_theta)
    random_radius = max_radius * numpy.sqrt(random.random())

    return numpy.array([
        random_radius * numpy.cos(random_theta),
        random_radius * numpy.sin(random_theta)
    ])

a = 2
b = 1

points = numpy.array([random_point(a, b) for _ in range(2000)])

plt.scatter(points[:,0], points[:,1])
plt.show()

Answer 4

我知道这是一个老问题，但我认为现有的答案都不够好。

我正在寻找完全相同问题的解决方案，并在 Google 的指导下找到了这里，发现所有现有答案都不是我想要的，所以我完全自己实现了自己的解决方案，使用此处找到的信息：https:// en.wikipedia.org/wiki/Ellipse

那么椭圆上的任意点都必须满足该方程，如何在椭圆内部制作一个点？

只需用 0 到 1 之间的两个随机数缩放 a 和 b。

我将在这里发布我的代码，我只是想提供帮助。

import math
import matplotlib.pyplot as plt
import random
from matplotlib.patches import Ellipse

a = 4
b = a*math.tan(math.radians((random.random()+0.5)/2*45))

def random_point(a, b):
    d = math.radians(random.random()*360)
    return (a * math.cos(d) * random.random(), b * math.sin(d) * random.random())

points = [random_point(a, b) for i in range(360)]

x, y = zip(*points)

fig = plt.figure(frameon=False)
ax = fig.add_subplot(111)
ax.set_axis_off()
ax.add_patch(Ellipse((0, 0), 2*a, 2*b, edgecolor='k', fc='None', lw=2))
ax.scatter(x, y)
fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
plt.axis('scaled')
plt.box(False)
ax = plt.gca()
ax.set_xlim([-a, a])
ax.set_ylim([-b, b])
plt.set_cmap('rainbow')
plt.show()

Answer 5

我们保持角度均匀随机，并让它稍后缩放。

def uniform_ellipse(N, a1, b2):

    rng = np.random.default_rng()
    
    #This line prevents overcrowding around the center
    x = (rng.uniform(size = N) * a1 ** 2) ** .5
    
    #Random angles from -pi to pi
    t1 = rng.uniform(-np.pi, np.pi, size = N)

    #Major axis points
    x1 = (x * np.cos(t1))
    
    #Scaling for the minor axis
    y1 = (x * np.sin(t1) * b2 / a1)
    
    #Plotting
    fig, ax = plt.subplots()
    ax.set_aspect('equal')
    ax.scatter(x1, y1, s = 4, color = 'k')
    ax.set_xticks([])
    ax.set_yticks([])
    plt.show()
    return x1, y1

其余的只是我试图检查结果是否确实均匀随机的测试。

我应用了两种不同的测试来检查结果是否确实是随机的。首先，我交叉检查了每个切片的面积以及该切片具有适当角度的点数。

def ellipse_slice(bin_count, a1, b2):
    area = lambda t1: np.arctan((b2/a1) * np.tan(t1))
    cum_area = np.vectorize(area)
    t_array = np.linspace(0, 0.5 * np.pi, bin_count + 1)
    slice_area = np.zeros(4 * bin_count)
    quadrant = np.diff(cum_area(t_array))
    slice_area[:bin_count] = quadrant
    slice_area[bin_count:2*bin_count] = quadrant[::-1]
    slice_area[2*bin_count:3*bin_count] = quadrant
    slice_area[3* bin_count:4*bin_count] = quadrant[::-1]
    t_array = np.linspace(-np.pi, np.pi, 4 * bin_count + 1)
    return slice_area * (b2 * a1 / 2) * np.pi

def angles_hist(bin_count, N, a1, b2):
    
    #Showing that angles are indeed random at the start
    fig, ax = plt.subplots()
    rng = np.random.default_rng()
    x = (rng.uniform(size = N) * a1 ** 2) ** .5
    t1 = rng.uniform(-np.pi, np.pi, size = N)
    bins = np.linspace(-np.pi, np.pi, bin_count + 1)
    n, bins, patches = ax.hist(t1, bins = bins, color = 'g', alpha = .6, density=True)    
    plt.show()
    
    #Checking if the angles now favors the major axis
    fig, ax = plt.subplots()
    x1 = (x * np.cos(t1))
    y1 = (x * np.sin(t1) * b2 / a1)
    tt = np.arctan2(x1, y1)
    bins = np.linspace(-np.pi, np.pi, bin_count + 1)
    n, bins, patches = ax.hist(tt, bins = bins, align='mid', color = 'g', alpha = .6, density=True)

    #plotting the area of each slice
    s1 = ellipse_slice(int(round(bin_count/4)), a1, b2) 
    ax.plot(bins[1:], s1)
    plt.show()
    return x1, y1

然后，我计算了到短轴和长轴的 epsilon 距离内的点数。我们期望点数应该与每个轴的长度成正比。

#Make sure to comment out the plotting part of uniform_ellipse
def point_counter(N, a1, b2, epsilon):
    x1, y1 = uniform_ellipse(N, a1, b2)
    cntx = sum(np.where(abs(x1)<epsilon, 1, 0))
    cnty = sum(np.where(abs(y1)<epsilon, 1, 0))
    return cntx, cnty

对于

a1 = 3, b2 = 1, N = 5000

我有

看起来还不错。这是切片的面积与该切片中点数的关系。这种转变可能是因为垃圾箱。

正如我所建议的，由于轴的缩放比例不同，角度会发生偏移。

最后，对于

N = 1000000, epsilon = 0.01

，短轴周围有 4260 个点，长轴周围有 12717 个点。