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