Numpy：生成一个二维高斯和pdf作为一个数组

Question

我正在尝试生成一个[600 x 600]numpy数组，它包含10个类似高斯的数组的和（每个数组都有一个随机生成的中心）。在

我试着用高斯滤波器来求和，但我确实用了高斯滤波器。即使使用num_centers=10也很慢，我可能需要和多达20个高斯数。在

这里有一个类似的问题，但似乎没有一个好的或结论性的答案，我不知道如何应用它来解决我的问题。 Sum of Gaussians into fast Numpy?

这是我尝试过的。在

import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt


num_centers = 10               # number of Gaussians to sum
sigma = 100                    # std. dev. of each Gaussian
result = np.zeros((600, 600))


for _ in range(num_centers):

    # Pick a random coordinate within the array as the center
    center = np.random.uniform(result.shape).astype(int)

    # Make array with 1 at the center and 0 everywhere else
    temp = np.zeros_like(result)
    temp[center[0], center[1]] = 1

    # Apply filter
    gaussian = gaussian_filter(temp, sigma)

    # Add to result
    result += gaussian


# Result should look like a contour map with several hills
plt.imshow(result * 1000)        # scale up to see the coloring
plt.show()

Answer 1

您可以消除循环，而是在每个中心创建一个值为1的数组，然后将gaussian_filter应用于该数组一次。所有的步骤都可以矢量化。在

这里有一个例子。我把sigma变小，这样就更容易区分中心了，我把宽度增加到800（没有什么特别的原因：）。在

import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt


num_centers = 10
sigma = 25
size = (600, 800)

impulses = np.zeros(size)

# rows and cols are the row and column indices of the centers
# of the gaussian peaks.
np.random.seed(123456)
rows, cols = np.unravel_index(np.random.choice(impulses.size, replace=False,
                                               size=num_centers),
                              impulses.shape)
impulses[rows, cols] = 1
# or use this if you want duplicates to sum:
# np.add.at(impulses, (rows, cols), 1)

# Filter impulses to create the result.
result = gaussian_filter(impulses, sigma, mode='nearest')

plt.imshow(result)
plt.show()

情节如下：

您可以尝试使用mode参数gaussian_filter，看看哪种模式最适合您。在

Answer 2

我不确定如何以并行方式处理随机高斯数组的创建，因为这是在代码中花费最多时间的。（我用timeit来确定这一点）。这是意料之中的，因为gaussian_filter是一个计算密集型函数。在

然而，我确实看到了在高斯数组上使用np.sum()的性能有了轻微的提高。这是因为调用np.sum()一次比从循环中调用+=更有效。在

示例

import numpy as np
from scipy.ndimage import gaussian_filter
import matplotlib.pyplot as plt


num_centers = 10               # number of Gaussians to sum
sigma = 100                    # std. dev. of each Gaussian
holder = np.zeros((num_centers, 600, 600))


for _ in range(num_centers):

    # Pick a random coordinate within the array as the center
    center = np.random.uniform(result.shape).astype(int)

    # Make array with 1 at the center and 0 everywhere else
    temp = np.zeros((600, 600))
    temp[center[0], center[1]] = 1

    # Apply filter
    holder[_] = gaussian_filter(temp, sigma)

result = np.sum(holder, axis=0)

# Result should look like a contour map with several hills
plt.imshow(result * 1000)        # scale up to see the coloring
plt.show()