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

possible duplicate of [Creating Gaussian filter of required length in python](http://stackoverflow.com/questions/11209115/creating-gaussian-filter-of-required-length-in-python) and http://astrolitterbox.blogspot.co.uk/2012/04/creating-discrete-gaussian-kernel-with.html — YXD, Jun 19 '13 at 12:13
I am using the code mentioned in the blog. I set `N = 2 and sigma = 1` and use 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)) / 2) print g / np.sqrt(2 * np.pi)` But the result obtained here is different form the one obtained using fspecial in matlab — Khushboo, Jun 19 '13 at 12:30
How is it different? What do you expect and what do you get? — interjay, Jun 19 '13 at 12:47

ali_m · Accepted Answer · 2013-09-10T17:04:41.623

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

Thats exactly what I was looking for. And you really made it simple. Thanks :) — Khushboo, Jun 20 '13 at 07:32
@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. — Moondra, Apr 04 '18 at 03:39
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. — Mong H. Ng, May 12 '19 at 16:32

score 14 · Answer 2 · edited Jan 13 '15 at 18:33

14

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

edited Jan 13 '15 at 18:33

J0e3gan

8,740
10
53
80

answered Jan 13 '15 at 17:58

sparklearner

403
3
7

Is there a reason this was downvoted? It seems to be correct (and more concise) – Moondra Apr 04 '18 at 03:35
I use this function and it has the nearest values to Matlab output. Thank you. – CC7052 Nov 15 '21 at 10:38

score 10 · Answer 3 · edited Oct 08 '17 at 15:27

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

score 3 · Answer 4 · answered Dec 11 '13 at 06:34

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.

score 1 · Answer 5 · answered Oct 23 '19 at 14:46

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

score 0 · Answer 6 · answered Mar 06 '16 at 20:47

Hi I think the problem is that for a gaussian filter the normalization factor depends on how many dimensions you used. So the filter looks like this formula
What you miss is the square of the normalization factor! And need to renormalize the whole matrix because of computing accuracy! The code is attached here:

def gaussian_filter(shape =(5,5), sigma=1):
    x, y = [edge /2 for edge in shape]
    grid = np.array([[((i**2+j**2)/(2.0*sigma**2)) for i in xrange(-x, x+1)] for j in xrange(-y, y+1)])
    g_filter = np.exp(-grid)/(2*np.pi*sigma**2)
    g_filter /= np.sum(g_filter)
    return g_filter
print gaussian_filter()

The output without normalized to sum of 1:

[[ 0.00291502  0.01306423  0.02153928  0.01306423  0.00291502]
 [ 0.01306423  0.05854983  0.09653235  0.05854983  0.01306423]
 [ 0.02153928  0.09653235  0.15915494  0.09653235  0.02153928]
 [ 0.01306423  0.05854983  0.09653235  0.05854983  0.01306423]
 [ 0.00291502  0.01306423  0.02153928  0.01306423  0.00291502]]

The output divided by np.sum(g_filter):

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

weiyixie · Answer 7 · 2020-01-25T08:43:38.060

here is to provide an nd-gaussian window generator:

def gen_gaussian_kernel(shape, mean, var):
    coors = [range(shape[d]) for d in range(len(shape))]
    k = np.zeros(shape=shape)
    cartesian_product = [[]]
    for coor in coors:
        cartesian_product = [x + [y] for x in cartesian_product for y in coor]
    for c in cartesian_product:
        s = 0
        for cc, m in zip(c,mean):
            s += (cc - m)**2
        k[tuple(c)] = np.exp(-s/(2*var))
    return k

this function will give you an unnormalized gaussian windows with given shape, center, and variance. for instance: gen_gaussian_kernel(shape=(3,3,3),mean=(1,1,1),var=1.0) output->

[[[ 0.22313016  0.36787944  0.22313016]
  [ 0.36787944  0.60653066  0.36787944]
  [ 0.22313016  0.36787944  0.22313016]]

 [[ 0.36787944  0.60653066  0.36787944]
  [ 0.60653066  1.          0.60653066]
  [ 0.36787944  0.60653066  0.36787944]]

 [[ 0.22313016  0.36787944  0.22313016]
  [ 0.36787944  0.60653066  0.36787944]
  [ 0.22313016  0.36787944  0.22313016]]]

This doesn't seem to work if the shape is just 2D...aka in the case I want to generate a single guassian kernel. — Will.Evo, Oct 02 '19 at 18:56
@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) — weiyixie, Jan 25 '20 at 08:45

score 0 · Answer 8 · answered Oct 13 '20 at 06:01

Using Gaussian PDF and assuming space invariant blur

def gaussian_kernel(sigma, size):
    mu = np.floor([size / 2, size / 2])
    size = int(size)
    kernel = np.zeros((size, size))
    for i in range(size):
        for j in range(size):
            kernel[i, j] = np.exp(-(0.5/(sigma*sigma)) * (np.square(i-mu[0]) + 
            np.square(j-mu[0]))) / np.sqrt(2*math.pi*sigma*sigma)```

    kernel = kernel/np.sum(kernel)
    return kernel

How to obtain a gaussian filter in python

8 Answers8

Linked