如何在大型数组中更快地生成随机点？

Question

我试图在3d胡椒形状内生成随机点，然而，当arr_size很大时，生成这些点需要太长时间。

当arr_size = (265,490,286)时，在3d形状内生成1000个随机点只需很短的时间，而arr_size为(30,30,30)需要很长的时间。

from matplotlib import pyplot as plt
import numpy as np


def create_bin_pepper(arr_size, center):
    coords = np.ogrid[:arr_size[0], :arr_size[1], :arr_size[2]]
    c = 10
    a1 = np.random.randint(low=5,high=10)
    b1 = np.random.randint(low=7,high=10)
    a2 = np.random.randint(low=5,high=10)
    b2 = np.random.randint(low=7,high=10)
    a3 = np.random.randint(low=5,high=10)
    b3 = np.random.randint(low=7,high=10)
    ellipse1 = ((np.square(coords[0] - center[0]))/np.square(a1) + (np.square(coords[1]-center[1]))/np.square(b1)  + (np.square(coords[2]-center[2]))/np.square(c)  <= 1)
    ellipse2 = ((np.square(coords[0] - center[0]-3))/np.square(a2) + (np.square(coords[1]-center[1]-5))/np.square(b2)  + (np.square(coords[2]-center[2]))/np.square(c)  <= 1)
    ellipse3 = ((np.square(coords[0] - center[0]+3))/np.square(a3) + (np.square(coords[1]-center[1]-5))/np.square(b3)  + (np.square(coords[2]-center[2]))/np.square(c)  <= 1)
    pepper = ellipse1|ellipse2|ellipse3
    pepper2 = np.where(pepper==1,230,pepper)

    for im in range(0,1000):
        #r2=1
        centre_x1 = np.random.randint(low=center[0]-a1+4,high=center[0]+a1-4)#low=11,high=20
        centre_y1 = np.random.randint(low=center[1]-b1+4,high=center[1]+b1-4)#low=15,high=23
        centre_z1 = np.random.randint(low=center[2]-c+4,high=center[2]+c-4)#low=10,high=20

        centre_x2 = np.random.randint(low=center[0]-a2+4,high=center[0]+a2-4)#low=11,high=20
        centre_y2 = np.random.randint(low=center[1]-b2+4,high=center[1]+b2-4)#low=15,high=23
        centre_z2 = np.random.randint(low=center[2]-c+4,high=center[2]+c-4)#low=10,high=20

        centre_x3 = np.random.randint(low=center[0]-a3+4,high=center[0]+a3-4)#low=11,high=20
        centre_y3 = np.random.randint(low=center[1]-b3+4,high=center[1]+b3-4)#low=15,high=23
        centre_z3 = np.random.randint(low=center[2]-c+4,high=center[2]+c-4)#low=10,high=20

        inside_ellipse1 = ((np.square(coords[0] - centre_x1))/np.square(a1) + (np.square(coords[1]-centre_y1))/np.square(b1)  + (np.square(coords[2]-centre_z1))/np.square(c)  <= (1/((np.square(a1))*(np.square(b1))*(np.square(c)))))
        inside_ellipse2 = ((np.square(coords[0] - centre_x2-3))/np.square(a2) + (np.square(coords[1]-centre_y2-5))/np.square(b2)  + (np.square(coords[2]-centre_z2))/np.square(c)  <= (1/((np.square(a2))*(np.square(b2))*(np.square(c)))))
        inside_ellipse3 = ((np.square(coords[0] - centre_x3+3))/np.square(a3) + (np.square(coords[1]-centre_y3-5))/np.square(b3)  + (np.square(coords[2]-centre_z3))/np.square(c)  <= (1/((np.square(a3))*(np.square(b3))*(np.square(c)))))

    pepper2 = inside_ellipse1 | inside_ellipse2 | inside_ellipse3 | pepper2
    pepper3 = np.where((pepper2!=230)&(pepper2!=0),160,pepper2)
    return pepper3

arr_size = (265,490,286)
sphere_center1 = (133,216,40)

pepper = create_bin_pepper(arr_size,sphere_center1)
axis = pepper[:,:,40]
plt.imshow(axis,cmap='gray')#,interpolation='bicubic'
plt.show()

Answer 1

通过在与arr_size形状相同的数组上使用np.random.rand，并使用椭圆和条件< 1000 / (area of ellipses)对其进行掩蔽，可以生成接近1000个点

from matplotlib import pyplot as plt
import numpy as np

POINTS = 1000

def ellipse(coords, center, offset):
    a = np.random.randint(low=5, high=10)
    b = np.random.randint(low=7, high=10)
    c = 10

    xs, ys, zs = coords
    cx, cy, cz = center
    ox, oy, oz = offset
    return ((xs - cx - ox) / a)**2 + ((ys - cy - oy) / b)**2 + ((zs - cz - oz) / c)**2  <= 1


def create_bin_pepper(arr_size, center):
    x, y, z = arr_size
    coords = np.ogrid[:x, :y, :z]
    ellipses = [ellipse(coords, center, offset) for offset in ((0, 0, 0), (3, 5, 0), (-3, 5, 0))]
    ellipses = np.logical_or.reduce(ellipses)

    area = ellipses.sum()

    random_points = np.where(ellipses, np.random.rand(*arr_size) < POINTS / area, 0)
    return random_points

arr_size = (300, 300, 300)
sphere_center = (150, 150, 150)

pepper = create_bin_pepper(arr_size, sphere_center)
print(pepper.sum())

它应该接近您需要生成的点数

Answer 2

欢迎使用StackOverflow！

看起来您是在每次循环迭代时重新分配centre_*和inside_ellipse*变量，而根本不使用以前的值。我假设这是基准测试中留下的，但如果不是这样，您可以简单地删除for im in range(0,1000)循环，并且已经实现了1000倍的加速。