efficiently generate "shifted" gaussian kernel in python

Question

I have a (very large) number of data points, each consisting of an x and y coordinate and a sigma-uncertainty (sigma is the same in both x and y directions; all three variables are floats). For each data-point I want to generate a 2d array on a standard grid, with probabilities that the the actual value is in that location.

For instance if x=5.0, y=5.0, sigma=1.0, on a (0,0)->(9,9) grid, I expect to generate:

   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.01,  0.02,  0.01,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.01,  0.06,  0.1 ,  0.06,  0.01,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.02,  0.1 ,  0.16,  0.1 ,  0.02,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.01,  0.06,  0.1 ,  0.06,  0.01,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.01,  0.02,  0.01,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
   [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ]]

Above was generated by creating a numpy array with zeroes, and [5,5] = 1, and then applying ndimage.filters.gaussian_filter with a sigma of 1. I feel that I can deal with non-integer x and y by distributing over nearby integer values and get a good approximation.

It feels however to be extreme overkill to get my resulting array this way, since scipy will have to take all values into account, not just the 1 in location [5, 5], even though they are all 0. It only takes 300us for a 64x64 grid, but still, I would likt to know if there is no more efficient way to get a X*Y numpy array with a gaussian kernel with arbitrary x, y and sigma.

Answer 1

A reasonably fast approach is to note that the Gaussian is separable, so you can calculate the 1D gaussian for x and y and then take the outer product:

import numpy as np
import matplotlib.pyplot as plt

x0, y0, sigma = 5.5, 4.2, 1.4

x, y = np.arange(9), np.arange(9)

gx = np.exp(-(x-x0)**2/(2*sigma**2))
gy = np.exp(-(y-y0)**2/(2*sigma**2))
g = np.outer(gx, gy)
g /= np.sum(g)  # normalize, if you want that

plt.imshow(g, interpolation="nearest", origin="lower")
plt.show()

Answer 2

@tom10's outer product answer is probably the best for this particular case. If you want to make a kernal out of an arbitrary function in two (or more) dimensions, you may want to look at np.indices or np.meshgrid.

For example:

def gaussian(x, mu=0, sigma=1):
    n = np.prod(sigma)*np.sqrt(2*np.pi)**len(x)
    return np.exp(-0.5*(((x-mu)/sigma)**2).sum(0))/n

gaussian(np.indices((10,10)), mu=5, sigma=1)

Which yields:

array([[ 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.001, 0.002, 0.001, 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.003, 0.013, 0.022, 0.013, 0.003, 0.   , 0.   ],
       [ 0.   , 0.   , 0.001, 0.013, 0.059, 0.097, 0.059, 0.013, 0.001, 0.   ],
       [ 0.   , 0.   , 0.002, 0.022, 0.097, 0.159, 0.097, 0.022, 0.002, 0.   ],
       [ 0.   , 0.   , 0.001, 0.013, 0.059, 0.097, 0.059, 0.013, 0.001, 0.   ],
       [ 0.   , 0.   , 0.   , 0.003, 0.013, 0.022, 0.013, 0.003, 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.001, 0.002, 0.001, 0.   , 0.   , 0.   ],
       [ 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.   ]])

For more flexibility, you can use np.meshgrid to control what the scale and scope of your domain is:

kern = gaussian(np.meshgrid(np.linspace(-10, 5), np.linspace(-2, 2)))

For this, kern.shape will be (50, 50) because 50 is the default length of np.linspace, and meshgrid is defining the x and y axes by the arrays passed to it. An equivalent way of doing this is np.mgrid[-10:5:50j, -2:2:50j]