# python – How to calculate a Gaussian kernel matrix efficiently in numpy?

## python – How to calculate a Gaussian kernel matrix efficiently in numpy?

I myself used the accepted answer for my image processing, but I find it (and the other answers) too dependent on other modules. Therefore, here is my compact solution:

``````import numpy as np

def gkern(l=5, sig=1.):

creates gaussian kernel with side length `l` and a sigma of `sig`

ax = np.linspace(-(l - 1) / 2., (l - 1) / 2., l)
gauss = np.exp(-0.5 * np.square(ax) / np.square(sig))
kernel = np.outer(gauss, gauss)
return kernel / np.sum(kernel)
``````

Edit: Changed arange to linspace to handle even side lengths

Edit: Use separability for faster computation, thank you Yves Daoust.

Do you want to use the Gaussian kernel for e.g. image smoothing? If so, theres a function `gaussian_filter()` in scipy:

This should work – while its still not 100% accurate, it attempts to account for the probability mass within each cell of the grid. I think that using the probability density at the midpoint of each cell is slightly less accurate, especially for small kernels. See https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm for an example.

``````import numpy as np
import scipy.stats as st

def gkern(kernlen=21, nsig=3):
Returns a 2D Gaussian kernel.

x = np.linspace(-nsig, nsig, kernlen+1)
kern1d = np.diff(st.norm.cdf(x))
kern2d = np.outer(kern1d, kern1d)
return kern2d/kern2d.sum()
``````

Testing it on the example in Figure 3 from the link:

``````gkern(5, 2.5)*273
``````

gives

``````array([[ 1.0278445 ,  4.10018648,  6.49510362,  4.10018648,  1.0278445 ],
[ 4.10018648, 16.35610171, 25.90969361, 16.35610171,  4.10018648],
[ 6.49510362, 25.90969361, 41.0435344 , 25.90969361,  6.49510362],
[ 4.10018648, 16.35610171, 25.90969361, 16.35610171,  4.10018648],
[ 1.0278445 ,  4.10018648,  6.49510362,  4.10018648,  1.0278445 ]])
``````

The original (accepted) answer below accepted is wrong
The square root is unnecessary, and the definition of the interval is incorrect.

``````import numpy as np
import scipy.stats as st

def gkern(kernlen=21, nsig=3):
Returns a 2D Gaussian kernel array.

interval = (2*nsig+1.)/(kernlen)
x = np.linspace(-nsig-interval/2., nsig+interval/2., kernlen+1)
kern1d = np.diff(st.norm.cdf(x))
kernel_raw = np.sqrt(np.outer(kern1d, kern1d))
kernel = kernel_raw/kernel_raw.sum()
return kernel
``````

#### python – How to calculate a Gaussian kernel matrix efficiently in numpy?

Im trying to improve on FuzzyDucks answer here. I think this approach is shorter and easier to understand. Here Im using `signal.scipy.gaussian` to get the 2D gaussian kernel.

``````import numpy as np
from scipy import signal

def gkern(kernlen=21, std=3):
Returns a 2D Gaussian kernel array.
gkern1d = signal.gaussian(kernlen, std=std).reshape(kernlen, 1)
gkern2d = np.outer(gkern1d, gkern1d)
return gkern2d
``````

Plotting it using `matplotlib.pyplot`:

``````import matplotlib.pyplot as plt
plt.imshow(gkern(21), interpolation=none)
``````