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

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

通过仔细考虑极坐标形式的定义，也可以在不使用拒绝采样的情况下生成椭圆内的点。根据维基百科，椭圆的极坐标形式由下式给出

\n\n

\n\n

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

\n\n

\n\n

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

\n\n

\n\n

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

\n\n

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

\n\n

下面是一些实现该算法的 Python 代码示例：

\n\n

import numpy\nimport matplotlib.pyplot as plt\nimport random\n\n# Returns theta in [-pi/2, 3pi/2]\ndef generate_theta(a, b):\n    u = random.random() / 4.0\n    theta = numpy.arctan(b/a * numpy.tan(2*numpy.pi*u))\n\n    v = random.random()\n    if v < 0.25:\n        return theta\n    elif v < 0.5:\n        return numpy.pi - theta\n    elif v < 0.75:\n        return numpy.pi + theta\n    else:\n        return -theta\n\ndef radius(a, b, theta):\n    return a * b / numpy.sqrt((b*numpy.cos(theta))**2 + (a*numpy.sin(theta))**2)\n\ndef random_point(a, b):\n    random_theta = generate_theta(a, b)\n    max_radius = radius(a, b, random_theta)\n    random_radius = max_radius * numpy.sqrt(random.random())\n\n    return numpy.array([\n        random_radius * numpy.cos(random_theta),\n        random_radius * numpy.sin(random_theta)\n    ])\n\na = 2\nb = 1\n\npoints = numpy.array([random_point(a, b) for _ in range(2000)])\n\nplt.scatter(points[:,0], points[:,1])\nplt.show()\n

\n\n

\n