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:

Updated answer

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)

Gaussian

Leave a Reply

Your email address will not be published.