在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

基于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

根据点的随机性,可能很简单,只需在磁盘中制作一个点网格,然后将每个点移动一些小但随机的量.