6

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.

Improve this question
5
  • I'm not sure I understand. Generating the kernel is the problem, not assigning it. Perhaps I should have been more clear. I want to generate a say 64 by 64 kernel for a 2d Gaussian with mean (6.7, 18.3) and sigma 12.6 in both directions. The mean and sigma are in the same units as the kernel. Commented Mar 9, 2015 at 20:07
  • I think I didn't understand you question correctly, and still don't. All the background and the way you are now doing it is fine, but please say clearly exactly what you want. My guess is that you just want a 9x9 Gaussian kernel where x, y, and sigma are specified but not necessarily integer valued. Is that right? Commented Mar 9, 2015 at 20:26
  • That is correct. And I want the "most efficient" way to do so, since I'll be doing it thousands of times a second with different values. Commented Mar 9, 2015 at 20:30
  • 1
    A reasonable approach to this is to note that the gaussian kernel is separable, so calculate the 1D for x, and the 1D for y, using the usual exp(-(x-x0)***2/sigma) and then take the outer product. This should be much faster than what you're going now anyway. Commented Mar 9, 2015 at 20:50
  • Nice!!! Never realised that, thanks! Well, I guess with "most efficient" I was just hoping that there was a scipy function that did exactly what I needed and I wouldn't have to worry that it could be optimised further. I guess that by your remark you didn't quite answer the question I asked, but did answer the question I should have asked :) Commented Mar 9, 2015 at 21:04

2 Answers 2

Reset to default
7

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

enter image description here

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

4

@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]

gaussian kernel

Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.