Numpy: Generating A 2d Sum Of Gaussians Pdf As An Array

I'm trying to generate a [600 x 600] numpy array that contains the sum of 10 Gaussian-like arrays (each with a randomly-generated center). I've tried using a Gaussian filter to gen

Solution 1:

You can eliminate the loop, and instead create an array with the value 1 at each center and then apply gaussian_filteronce to this array. All the steps can be vectorized.

Here's an example. I made sigma smaller so it was easier to distinguish the centers, and I increased the width to 800 (for no particular reason :).

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

Here's the plot:

You can experiment with the mode argument of gaussian_filter to see which mode works best for you.

Solution 2:

Im not sure how you can handle the creation of the random gaussian arrays in a parrallel fashion, as that is what is taking the most time in your code. (I used timeit to determine this). This is to be expected, as gaussian_filter is a computationally intensive function.

However, I did see a slight performance boost by using np.sum() on an array of gaussians. This is because calling np.sum() once is more efficient than calling += from within a loop.

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