Numpy：生成一个二维Gaussians 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()一次比在循环中调用+=更有效。

Example

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()