How to obtain a gaussian filter in python

Question

I am using python to create a gaussian filter of size 5x5. I saw this post here where they talk about a similar thing but I didn't find the exact way to get equivalent python code to matlab function fspecial('gaussian', f_wid, sigma) Is there any other way to do it? I tried using the following code :

size = 2
sizey = None
size = int(size)
if not sizey:
    sizey = size
else:
    sizey = int(sizey)
x, y = scipy.mgrid[-size: size + 1, -sizey: sizey + 1]
g = scipy.exp(- (x ** 2/float(size) + y ** 2 / float(sizey)))
print g / np.sqrt(2 * np.pi)

The output obtained is

[[ 0.00730688  0.03274718  0.05399097  0.03274718  0.00730688]
 [ 0.03274718  0.14676266  0.24197072  0.14676266  0.03274718]
 [ 0.05399097  0.24197072  0.39894228  0.24197072  0.05399097]
 [ 0.03274718  0.14676266  0.24197072  0.14676266  0.03274718]
 [ 0.00730688  0.03274718  0.05399097  0.03274718  0.00730688]]

What I want is something like this:

   0.0029690   0.0133062   0.0219382   0.0133062   0.0029690
   0.0133062   0.0596343   0.0983203   0.0596343   0.0133062
   0.0219382   0.0983203   0.1621028   0.0983203   0.0219382
   0.0133062   0.0596343   0.0983203   0.0596343   0.0133062
   0.0029690   0.0133062   0.0219382   0.0133062   0.0029690

Solution 1

In general terms if you really care about getting the the exact same result as MATLAB, the easiest way to achieve this is often by looking directly at the source of the MATLAB function.

In this case, edit fspecial:

...
  case 'gaussian' % Gaussian filter

     siz   = (p2-1)/2;
     std   = p3;

     [x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1));
     arg   = -(x.*x + y.*y)/(2*std*std);

     h     = exp(arg);
     h(h<eps*max(h(:))) = 0;

     sumh = sum(h(:));
     if sumh ~= 0,
       h  = h/sumh;
     end;
...

Pretty simple, eh? It's <10mins work to port this to Python:

import numpy as np

def matlab_style_gauss2D(shape=(3,3),sigma=0.5):
    """
    2D gaussian mask - should give the same result as MATLAB's
    fspecial('gaussian',[shape],[sigma])
    """
    m,n = [(ss-1.)/2. for ss in shape]
    y,x = np.ogrid[-m:m+1,-n:n+1]
    h = np.exp( -(x*x + y*y) / (2.*sigma*sigma) )
    h[ h < np.finfo(h.dtype).eps*h.max() ] = 0
    sumh = h.sum()
    if sumh != 0:
        h /= sumh
    return h

This gives me the same answer as fspecial to within rounding error:

 >> fspecial('gaussian',5,1)

 0.002969     0.013306     0.021938     0.013306     0.002969
 0.013306     0.059634      0.09832     0.059634     0.013306
 0.021938      0.09832       0.1621      0.09832     0.021938
 0.013306     0.059634      0.09832     0.059634     0.013306
 0.002969     0.013306     0.021938     0.013306     0.002969

 : matlab_style_gauss2D((5,5),1)

array([[ 0.002969,  0.013306,  0.021938,  0.013306,  0.002969],
       [ 0.013306,  0.059634,  0.09832 ,  0.059634,  0.013306],
       [ 0.021938,  0.09832 ,  0.162103,  0.09832 ,  0.021938],
       [ 0.013306,  0.059634,  0.09832 ,  0.059634,  0.013306],
       [ 0.002969,  0.013306,  0.021938,  0.013306,  0.002969]])

Solution 2

I found similar solution for this problem:

def fspecial_gauss(size, sigma):

    """Function to mimic the 'fspecial' gaussian MATLAB function
    """

    x, y = numpy.mgrid[-size//2 + 1:size//2 + 1, -size//2 + 1:size//2 + 1]
    g = numpy.exp(-((x**2 + y**2)/(2.0*sigma**2)))
    return g/g.sum()

Solution 3

You could try this too (as product of 2 independent 1D Gaussian random variables) to obtain a 2D Gaussian Kernel:

from numpy import pi, exp, sqrt
s, k = 1, 2 #  generate a (2k+1)x(2k+1) gaussian kernel with mean=0 and sigma = s
probs = [exp(-z*z/(2*s*s))/sqrt(2*pi*s*s) for z in range(-k,k+1)] 
kernel = np.outer(probs, probs)
print kernel

#[[ 0.00291502  0.00792386  0.02153928  0.00792386  0.00291502]
#[ 0.00792386  0.02153928  0.05854983  0.02153928  0.00792386]
#[ 0.02153928  0.05854983  0.15915494  0.05854983  0.02153928]
#[ 0.00792386  0.02153928  0.05854983  0.02153928  0.00792386]
#[ 0.00291502  0.00792386  0.02153928  0.00792386  0.00291502]]

import matplotlib.pylab as plt
plt.imshow(kernel)
plt.colorbar()
plt.show()

Solution 4

This function implements functionality similar to fspecial in matlab

http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.get_window.html from scipy import signal

>>>signal.get_window(('gaussian',2),3)
>>>array([ 0.8824969,  1.       ,  0.8824969])

This function appears to generate only 1D kernels

I guess you could implement code to generate a Gaussian mask yourself as well as other have pointed out.

Solution 5

Hey, I think this might help you

import numpy as np
import cv2

def gaussian_kernel(dimension_x, dimension_y, sigma_x, sigma_y):
    x = cv2.getGaussianKernel(dimension_x, sigma_x)
    y = cv2.getGaussianKernel(dimension_y, sigma_y)
    kernel = x.dot(y.T)
    return kernel
g_kernel = gaussian_kernel(5, 5, 1, 1)
print(g_kernel)

[[0.00296902 0.01330621 0.02193823 0.01330621 0.00296902]
 [0.01330621 0.0596343  0.09832033 0.0596343  0.01330621]
 [0.02193823 0.09832033 0.16210282 0.09832033 0.02193823]
 [0.01330621 0.0596343  0.09832033 0.0596343  0.01330621]
 [0.00296902 0.01330621 0.02193823 0.01330621 0.00296902]]

View more solutions

Comments

[Khushboo]

Thats exactly what I was looking for. And you really made it simple. Thanks :)

[Moondra]

Is there a reason this was downvoted? It seems to be correct (and more concise)

[Moondra]

@ali_m I don't quite understand this line h[ h < np.finfo(h.dtype).eps*h.max() ] = 0 I took a look at the docs and it seems to prevent machine limit errors(?). Is the d.type argument necessary? And why do you use .eps and multiply it with h.max. Sorry to bring up such an old post.

[Mong H. Ng]

It is also possible to use cv2 and the following 2 liners: u = c.getGaussianKernel(5, 1); kernel = np.outer(u, u). I got the same results as above.

[Will.Evo]

This doesn't seem to work if the shape is just 2D...aka in the case I want to generate a single guassian kernel.

[weiyixie]

@Will.Evo it works for 2d as well? I don't see why it is not working for 2d case. try run gen_gaussian_kernel((5, 5), (1.0, 1.0), 1.0)

[CC7052]

I use this function and it has the nearest values to Matlab output. Thank you.