用python中的点均匀随机填充磁盘的方法

Question

我有一个应用程序需要以准随机方式填充“n”个点的磁盘。我希望这些点有点随机，但在磁盘上仍然或多或少有规律的密度。

我目前的方法是放置一个点，检查它是否在磁盘内，然后检查它是否也离所有其他已保存的点足够远。我的代码如下：

import os
import random
import math

# ------------------------------------------------ #
# geometric constants
center_x = -1188.2
center_y = -576.9
center_z = -3638.3

disk_distance = 2.0*5465.6
disk_diam = 5465.6

# ------------------------------------------------ #

pts_per_disk = 256
closeness_criteria = 200.0
min_closeness_criteria = disk_diam/closeness_criteria

disk_center = [(center_x-disk_distance),center_y,center_z]
pts_in_disk = []
while len(pts_in_disk) < (pts_per_disk):
    potential_pt_x = disk_center[0]
    potential_pt_dy = random.uniform(-disk_diam/2.0, disk_diam/2.0)
    potential_pt_y = disk_center[1]+potential_pt_dy
    potential_pt_dz = random.uniform(-disk_diam/2.0, disk_diam/2.0)
    potential_pt_z = disk_center[2]+potential_pt_dz
    potential_pt_rad = math.sqrt((potential_pt_dy)**2+(potential_pt_dz)**2)

    if potential_pt_rad < (disk_diam/2.0):
        far_enough_away = True
        for pt in pts_in_disk:
            if math.sqrt((potential_pt_x - pt[0])**2+(potential_pt_y - pt[1])**2+(potential_pt_z - pt[2])**2) > min_closeness_criteria:
                pass
            else:
                far_enough_away = False
                break
        if far_enough_away:
            pts_in_disk.append([potential_pt_x,potential_pt_y,potential_pt_z])

outfile_name = "pt_locs_x_lo_"+str(pts_per_disk)+"_pts.txt"
outfile = open(outfile_name,'w')
for pt in pts_in_disk:
    outfile.write(" ".join([("%.5f" % (pt[0]/1000.0)),("%.5f" % (pt[1]/1000.0)),("%.5f" % (pt[2]/1000.0))])+'\n')
outfile.close()

为了获得最均匀的点密度，我所做的基本上是使用另一个脚本迭代地运行此脚本，每次连续迭代都会降低“接近度”标准。在某些时候，脚本无法完成，我只是使用上次成功迭代的点。

所以我的问题相当广泛：有没有更好的方法来做到这一点？我的方法目前还可以，但我的直觉告诉我有更好的方法来生成这样的点场。

输出的图示如下图，一个具有高接近度标准，另一个具有“最低找到”接近度标准（我想要的）。

Answer 1

基于Disk Point Picking from MathWorld的简单解决方案：

import numpy as np
import matplotlib.pyplot as plt

n = 1000
r = np.random.uniform(low=0, high=1, size=n)  # radius
theta = np.random.uniform(low=0, high=2*np.pi, size=n)  # angle

x = np.sqrt(r) * np.cos(theta)
y = np.sqrt(r) * np.sin(theta)

# for plotting circle line:
a = np.linspace(0, 2*np.pi, 500)
cx,cy = np.cos(a), np.sin(a)

fg, ax = plt.subplots(1, 1)
ax.plot(cx, cy,'-', alpha=.5) # draw unit circle line
ax.plot(x, y, '.') # plot random points
ax.axis('equal')
ax.grid(True)
fg.canvas.draw()
plt.show()

它给。

或者，您也可以创建一个规则网格并随机扭曲它：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri


n = 20
tt = np.linspace(-1, 1, n)
xx, yy = np.meshgrid(tt, tt)  # create unit square grid
s_x, s_y  = xx.ravel(), yy.ravel()
ii = np.argwhere(s_x**2 + s_y**2 <= 1).ravel() # mask off unwanted points
x, y = s_x[ii], s_y[ii]
triang = tri.Triangulation(x, y) # create triangluar grid


# distort the grid
g = .5 # distortion factor
rx = x + np.random.uniform(low=-g/n, high=g/n, size=x.shape)
ry = y + np.random.uniform(low=-g/n, high=g/n, size=y.shape)

rtri = tri.Triangulation(rx, ry, triang.triangles)  # distorted grid

# for circle:
a = np.linspace(0, 2*np.pi, 500)
cx,cy = np.cos(a), np.sin(a)


fg, ax = plt.subplots(1, 1)
ax.plot(cx, cy,'k-', alpha=.2) # circle line
ax.triplot(triang, "g-", alpha=.4)
ax.triplot(rtri, 'b-', alpha=.5)
ax.axis('equal')
ax.grid(True)
fg.canvas.draw()
plt.show()

它给

三角形只是为了可视化。明显的缺点是，根据您对网格的选择，无论是在中间还是在边界上（如此处所示），由于网格离散化，都会或多或少出现较大的“洞”。

Answer 2

如果您有一个像圆盘（圆）这样的定义区域，您希望在其中生成随机点，最好使用圆的方程并限制半径：

x^2 + y^2 = r^2  (0 < r < R)

或参数化为两个变量

cos(a) = x/r
sin(a) = y/r
sin^2(a) + cos^2(a) = 1

要生成类似低密度的伪随机分布，您应该采用以下方法：

对于r 和a 的随机分布范围，选择n 个点。

这使您可以生成大致符合您的密度标准的分布。

要理解为什么会这样，想象你的圆圈首先被分成长度为dr 的小环，现在想象你的圆圈被分成角度为da 的饼片。您的随机性现在在圆圈周围的整个盒装区域上具有相同的概率。如果您在整个圈子中划分允许的随机区域，您将在整个圈子周围获得更均匀的分布，并且各个区域的随机变化较小，从而为您提供您所追求的伪随机外观。

现在您的工作就是为每个给定区域生成n 点。您将希望n 依赖于r，因为当您移出圆圈时，每个分区的区域都会发生变化。您可以将其与每个空间带来的确切面积变化成比例：

对于第 n 到 n+1 环：

d(Area,n,n-1) = Area(n) - Area(n-1)

任何给定环的面积是：

Area = pi*(dr*n)^2 - pi*(dr*(n-1))

所以区别就变成了：

d(Area,n,n-1) = [pi*(dr*n)^2 - pi*(dr*(n-1))^2] - [pi*(dr*(n-1))^2 - pi*(dr*(n-2))^2]
d(Area,n,n-1) = pi*[(dr*n)^2 - 2*(dr*(n-1))^2 + (dr*(n-2))^2]

您可以对此进行阐述，以深入了解n 应该增加多少，但仅猜测某个百分比增加 (30%) 之类的可能会更快。

我提供的示例是一个小子集，减少 da 和 dr 将显着改善您的结果。

这里是一些用于生成这些点的粗略代码：

import random
import math

R = 10.
n_rings = 10.
n_angles = 10.
dr = 10./n_rings
da = 2*math.pi/n_angles

base_points_per_division = 3
increase_per_level = 1.1
points = []
ring = 0
while ring < n_rings:
    angle = 0
    while angle < n_angles:
        for i in xrange(int(base_points_per_division)):
             ra = angle*da + da*math.random()
             rr = r*dr + dr*random.random()
             x = rr*math.cos(ra)
             y = rr*math.sin(ra)
             points.append((x,y))
        angle += 1
    base_points_per_division = base_points_per_division*increase_per_level
    ring += 1

我用参数测试了一下：

n_rings = 20
n_angles = 20
base_points = .9
increase_per_level = 1.1

并得到以下结果：

它看起来比您提供的图像更密集，但我想进一步调整这些变量可能是有益的。

您可以添加一个额外的部分，通过计算每个环的点数来适当地缩放密度。

points_per_ring = 密度math.pi(dr**2)*(2*n+1) points_per_division = points_per_ring/n_angles

这将提供一个更好的规模分布。

density = .03
points = []
ring = 0
while ring < n_rings:
    angle = 0
    base_points_per_division = density*math.pi*(dr**2)*(2*ring+1)/n_angles
    while angle < n_angles:
        for i in xrange(int(base_points_per_division)):
             ra = angle*da + min(da,da*random.random())
             rr = ring*dr + dr*random.random()
             x = rr*math.cos(ra)
             y = rr*math.sin(ra)
             points.append((x,y))
        angle += 1
    ring += 1

使用以下参数给出更好的结果

R = 1.
n_rings = 10.
n_angles = 10.
density = 10/(dr*da)   # ~ ten points per unit area

用图表...

为了好玩，您可以绘制分区图，查看它与您的分布匹配的程度并进行调整。

Answer 3

根据点需要的随机性，只需在磁盘内制作一个点网格，然后将每个点替换一些很小但随机的量，这可能很简单。

Answer 4

您可能想要更多的随机性，但如果您只是想在圆盘上填充看起来均匀且不在明显网格上的点分布，您可以尝试使用随机相位的螺旋。

import math
import random
import pylab

n = 300

alpha = math.pi * (3 - math.sqrt(5))    # the "golden angle"
phase = random.random() * 2 * math.pi
points = []
for k in xrange(n):
  theta = k * alpha + phase
  r = math.sqrt(float(k)/n)
  points.append((r * math.cos(theta), r * math.sin(theta)))

pylab.scatter(*zip(*points))
pylab.show()

Answer 5

概率论确保rejection method 是一种合适的方法 在圆盘内生成均匀分布的点 D(0,r)，以原点为中心，半径为 r。即在正方形 [-r,r] x [-r,r] 内生成点，直到一个点落入圆盘内：

do{
  generate P in [-r,r]x[-r,r];
  }while(P[0]**2+P[1]**2>r);
  return P;

unif_rnd_disk 是实现此拒绝方法的生成器函数：

import matplotlib.pyplot as plt
import numpy as np
import itertools

def unif_rnd_disk(r=1.0):
   pt=np.zeros(2) 
   while True:
      yield pt  
      while True: 
        pt=-r+2*r*np.random.random(2)  
        if (pt[0]**2+pt[1]**2<=r):
            break


G=unif_rnd_disk()# generator of points in disk D(0,r=1)
X,Y=zip(*[pt for pt in itertools.islice(G, 1, 1000)])

plt.scatter(X, Y, color='r', s=3)
plt.axis('equal')

如果我们想在以 C(a,b) 为中心的圆盘中生成点，我们必须对圆盘 D(0,r) 中的点应用平移：

C=[2.0, -3.5]
plt.scatter(C[0]+np.array(X), C[1]+np.array(Y), color='r', s=3)
plt.axis('equal')