How to calculate a Gaussian kernel matrix efficiently in numpy?

Question

def GaussianMatrix(X,sigma):
    row,col=X.shape
    GassMatrix=np.zeros(shape=(row,row))
    X=np.asarray(X)
    i=0
    for v_i in X:
        j=0
        for v_j in X:
            GassMatrix[i,j]=Gaussian(v_i.T,v_j.T,sigma)
            j+=1
        i+=1
    return GassMatrix
def Gaussian(x,z,sigma):
    return np.exp((-(np.linalg.norm(x-z)**2))/(2*sigma**2))

This is my current way. Is there any way I can use matrix operation to do this? X is the data points.

Answer 1

Combining multiple suggestions results in sampling a Gaussian kernel on a given rectangular window size. Sigma is determined per axis of the window. Obviously, the following could be formulated more densely, however for the sake of readability and clearness:

#!/usr/bin/env python3
import numpy as np
import cv2
import scipy.signal

import matplotlib.pyplot as plt

def GaussKernel(windowShape):
    """
    https://stackoverflow.com/a/43346070/1150303
    https://dsp.stackexchange.com/a/74157
    """
    l = np.array(windowShape[:2])
    sig = (l - 1) / 4.
    
    x = np.linspace(-(l[0] - 1) / 2., (l[0] - 1) / 2., l[0])
    y = np.linspace(-(l[1] - 1) / 2., (l[1] - 1) / 2., l[1])    
    
    gaussX = np.exp(-0.5 * np.square(x) / np.square(sig[0]))
    gaussY = np.exp(-0.5 * np.square(y) / np.square(sig[1]))    
    kernel = np.outer(gaussX, gaussY)
    return kernel / np.sum(kernel)
    

for rows in range(31, 122, 30):
    for cols in range(31, 122, 30):
        kernel = GaussKernel( (rows, cols) )
        vis = (kernel / kernel.max() * 255).astype(np.uint8)
        vis = cv2.applyColorMap(vis, cv2.COLORMAP_JET)
        cv2.imshow(f"kernel rows {rows}, cols {cols}", vis)
cv2.waitKey(0)
        

Answer 2

Adapting the accepted answer by FuzzyDuck to match the results of this website: http://dev.theomader.com/gaussian-kernel-calculator/ I now present this definition to you:

import numpy as np
import scipy.stats as st

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

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

print(gkern(kernlen=5, sig=1))

output:

[[0.003765   0.015019   0.02379159 0.015019   0.003765  ]
 [0.015019   0.05991246 0.0949073  0.05991246 0.015019  ]
 [0.02379159 0.0949073  0.15034262 0.0949073  0.02379159]
 [0.015019   0.05991246 0.0949073  0.05991246 0.015019  ]
 [0.003765   0.015019   0.02379159 0.015019   0.003765  ]]

Answer 3

If you are a computer vision engineer and you need heatmap for a particular point as Gaussian distribution (especially for keypoint detection on image)

def gaussian_heatmap(center = (2, 2), image_size = (10, 10), sig = 1):
    """
    It produces single gaussian at expected center
    :param center:  the mean position (X, Y) - where high value expected
    :param image_size: The total image size (width, height)
    :param sig: The sigma value
    :return:
    """
    x_axis = np.linspace(0, image_size[0]-1, image_size[0]) - center[0]
    y_axis = np.linspace(0, image_size[1]-1, image_size[1]) - center[1]
    xx, yy = np.meshgrid(x_axis, y_axis)
    kernel = np.exp(-0.5 * (np.square(xx) + np.square(yy)) / np.square(sig))
    return kernel

The usage and output

kernel = gaussian_heatmap(center = (2, 2), image_size = (10, 10), sig = 1)
plt.imshow(kernel)
print("max at :", np.unravel_index(kernel.argmax(), kernel.shape))
print("kernel shape", kernel.shape)

max at : (2, 2)

kernel shape (10, 10)

kernel = gaussian_heatmap(center = (25, 40), image_size = (100, 50), sig = 5)
plt.imshow(kernel)
print("max at :", np.unravel_index(kernel.argmax(), kernel.shape))
print("kernel shape", kernel.shape)

max at : (40, 25)

kernel shape (50, 100)

Answer 4

I tried using numpy only. Here is the code

def get_gauss_kernel(size=3,sigma=1):
    center=(int)(size/2)
    kernel=np.zeros((size,size))
    for i in range(size):
       for j in range(size):
          diff_sq = (i-center)**2+(j-center)**2
          kernel[i,j]=np.exp(-diff_sq/(2*sigma**2))
    return kernel/np.sum(kernel)

You can visualise the result using:

plt.imshow(get_gauss_kernel(5,1))

Answer 5

I took a similar approach to Nils Werner's answer -- since convolution of any kernel with a Kronecker delta results in the kernel itself centered around that Kronecker delta -- but I made it slightly more general to deal with both odd and even dimensions. In three lines:

import scipy.ndimage as scim

def gaussian_kernel(dimension: int, sigma: float):
    dirac = np.zeros((dimension,dimension))
    dirac[(dimension-1)//2:dimension//2+1, (dimension-1)//2:dimension//2+1] = 1.0 / (1 + 3 * ((dimension + 1) % 2))
    return scim.gaussian_filter(dirac, sigma=sigma)

The second line creates either a single 1.0 in the middle of the matrix (if the dimension is odd), or a square of four 0.25 elements (if the dimension is even). The division could be moved to the third line too; the result is normalised either way.

For those who like to have the kernel the matrix with one (odd) or four (even) 1.0 element(s) in the middle instead of normalisation, this works:

import scipy.ndimage as scim

def gaussian_kernel(dimension: int, sigma: float, ones_in_the_middle=False):
    dirac = np.zeros((dimension,dimension))
    dirac[(dimension-1)//2:dimension//2+1, (dimension-1)//2:dimension//2+1] = 1.0
    kernel = scim.gaussian_filter(dirac, sigma=sigma)
    divisor = kernel[(dimension-1)//2, (dimension-1)//2] if ones_in_the_middle else 1 + 3 * ((dimension + 1) % 2)
    return kernel/divisor

Answer 6

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.

Answer 7

Yet another implementation.

This is normalized so that for sigma > 1 and sufficiently large win_size, the total sum of the kernel elements equals 1.

def gaussian_kernel(win_size, sigma):
    t = np.arange(win_size)
    x, y = np.meshgrid(t, t)
    o = (win_size - 1) / 2
    r = np.sqrt((x - o)**2 + (y - o)**2)
    scale = 1 / (sigma**2 * 2 * np.pi)
    return scale * np.exp(-0.5 * (r / sigma)**2)

To generate a 5x5 kernel:

gaussian_kernel(win_size=5, sigma=1)

Answer 8

A good way to do that is to use the gaussian_filter function to recover the kernel. For instance:

indicatrice = np.zeros((5,5))
indicatrice[2,2] = 1
gaussian_kernel = gaussian_filter(indicatrice, sigma=1)
gaussian_kernel/=gaussian_kernel[2,2]

This gives

array[[0.02144593, 0.08887207, 0.14644428, 0.08887207, 0.02144593],
       [0.08887207, 0.36828649, 0.60686612, 0.36828649, 0.08887207],
       [0.14644428, 0.60686612, 1.        , 0.60686612, 0.14644428],
       [0.08887207, 0.36828649, 0.60686612, 0.36828649, 0.08887207],
       [0.02144593, 0.08887207, 0.14644428, 0.08887207, 0.02144593]]

Answer 9

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

Updated answer

This should work - while it's 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

Answer 10

As I didn't find what I was looking for, I coded my own one-liner. You can modify it accordingly (according to the dimensions and the standard deviation).

Here is the one-liner function for a 3x5 patch for example.

from scipy import signal

def gaussian2D(patchHeight, patchWidth, stdHeight=1, stdWidth=1):
    gaussianWindow = signal.gaussian(patchHeight, stdHeight).reshape(-1, 1)@signal.gaussian(patchWidth, stdWidth).reshape(1, -1)
    return gaussianWindow

print(gaussian2D(3, 5))

You get an output like this:

[[0.082085   0.36787944 0.60653066 0.36787944 0.082085  ]
[0.13533528  0.60653066 1.         0.60653066 0.13533528]
[0.082085   0.36787944 0.60653066 0.36787944 0.082085  ]]

You can read more about scipy's Gaussian here.

Answer 11

A 2D gaussian kernel matrix can be computed with numpy broadcasting,

def gaussian_kernel(size=21, sigma=3):
    """Returns a 2D Gaussian kernel.
    Parameters
    ----------
    size : float, the kernel size (will be square)

    sigma : float, the sigma Gaussian parameter

    Returns
    -------
    out : array, shape = (size, size)
      an array with the centered gaussian kernel
    """
    x = np.linspace(- (size // 2), size // 2)
    x /= np.sqrt(2)*sigma
    x2 = x**2
    kernel = np.exp(- x2[:, None] - x2[None, :])
    return kernel / kernel.sum()

For small kernel sizes this should be reasonably fast.

Note: this makes changing the sigma parameter easier with respect to the accepted answer.

Answer 12

Building up on Teddy Hartanto's answer. You can just calculate your own one dimensional Gaussian functions and then use np.outer to calculate the two dimensional one. Very fast and efficient way.

With the code below you can also use different Sigmas for every dimension

import numpy as np
def generate_gaussian_mask(shape, sigma, sigma_y=None):
    if sigma_y==None:
        sigma_y=sigma
    rows, cols = shape

    def get_gaussian_fct(size, sigma):
        fct_gaus_x = np.linspace(0,size,size)
        fct_gaus_x = fct_gaus_x-size/2
        fct_gaus_x = fct_gaus_x**2
        fct_gaus_x = fct_gaus_x/(2*sigma**2)
        fct_gaus_x = np.exp(-fct_gaus_x)
        return fct_gaus_x

    mask = np.outer(get_gaussian_fct(rows,sigma), get_gaussian_fct(cols,sigma_y))
    return mask

Answer 13

I'm trying to improve on FuzzyDuck's answer here. I think this approach is shorter and easier to understand. Here I'm 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')

Answer 14

You may simply gaussian-filter a simple 2D dirac function, the result is then the filter function that was being used:

import numpy as np
import scipy.ndimage.filters as fi

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

    # create nxn zeros
    inp = np.zeros((kernlen, kernlen))
    # set element at the middle to one, a dirac delta
    inp[kernlen//2, kernlen//2] = 1
    # gaussian-smooth the dirac, resulting in a gaussian filter mask
    return fi.gaussian_filter(inp, nsig)

Answer 15

linalg.norm takes an axis parameter. With a little experimentation I found I could calculate the norm for all combinations of rows with

np.linalg.norm(x[None,:,:]-x[:,None,:],axis=2)

It expands x into a 3d array of all differences, and takes the norm on the last dimension.

So I can apply this to your code by adding the axis parameter to your Gaussian:

def Gaussian(x,z,sigma,axis=None):
    return np.exp((-(np.linalg.norm(x-z, axis=axis)**2))/(2*sigma**2))

x=np.arange(12).reshape(3,4)
GaussianMatrix(x,1)

produces

array([[  1.00000000e+00,   1.26641655e-14,   2.57220937e-56],
       [  1.26641655e-14,   1.00000000e+00,   1.26641655e-14],
       [  2.57220937e-56,   1.26641655e-14,   1.00000000e+00]])

Matching:

Gaussian(x[None,:,:],x[:,None,:],1,axis=2)

array([[  1.00000000e+00,   1.26641655e-14,   2.57220937e-56],
       [  1.26641655e-14,   1.00000000e+00,   1.26641655e-14],
       [  2.57220937e-56,   1.26641655e-14,   1.00000000e+00]])