CODEWITHSUNDEEP
pythonarraysnumpyvectorizationgaussian
Kanad Khanna
Kanad Khanna

Reputation: 33

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 generate the individual Gaussian-like arrays, then summing them up, but I'm sure there's a vectorized way to approach this. Even with num_centers=10 it's slow, and I might need to sum as many as 20 Gaussians.

There is a similar question here, but it doesn't seem to have a good or conclusive answer and I'm not sure how to apply it to my problem. Sum of Gaussians into fast Numpy?

Here's what I have tried.

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

Upvotes: 3

Views: 868

Answers (2)

Warren Weckesser
Warren Weckesser

Reputation: 114976

You can eliminate the loop, and instead create an array with the value 1 at each center and then apply gaussian_filter once 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:

plot

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

Upvotes: 2

nickanna_42
nickanna_42

Reputation: 179

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

Upvotes: 0

Related Questions