How to obtain a gaussian filter in python

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 npdef 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]])

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, sqrts, k = 1, 2 #  generate a (2k+1)x(2k+1) gaussian kernel with mean=0 and sigma = sprobs = [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 pltplt.imshow(kernel)plt.colorbar()plt.show()

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