Как эффективно вычислить гауссовскую матрицу ядра в numpy?

Ответ 1

Вы хотите использовать ядро Гаусса, например, для сглаживания изображений? Если это так, в scipy есть функция gaussian_filter():

Обновленный ответ

Это должно сработать - хотя оно все еще не является точным на 100%, оно пытается учесть массу вероятности в каждой ячейке сетки. Я думаю, что использование плотности вероятности в средней точке каждой ячейки немного менее точно, особенно для небольших ядер. См. Пример https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm.

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

Тестируем его на примере на рисунке 3 по ссылке:

gkern(5, 2.5)*273

дает

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

Исходный (принятый) ответ ниже принятый является неправильным Квадратный корень не нужен, и определение интервала является неправильным.

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

Ответ 2

Вы можете просто gaussian-фильтровать простую функцию 2D dirac, результатом будет функция фильтра, которая используется:

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)

Ответ 3

Я сам использовал принятый ответ для обработки изображений, но нахожу его (и другие ответы) слишком зависимым от других модулей. Поэтому вот мое компактное решение:

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)
    xx, yy = np.meshgrid(ax, ax)

    kernel = np.exp(-0.5 * (np.square(xx) + np.square(yy)) / np.square(sig))

    return kernel / np.sum(kernel)

Редактировать: изменено значение arange на linspace для обработки равных сторон

Ответ 4

Я пытаюсь улучшить ответ FuzzyDuck здесь. Я думаю, что этот подход короче и проще для понимания. Здесь я использую signal.scipy.gaussian для получения двумерного гауссовского ядра.

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

matplotlib.pyplot графика с использованием matplotlib.pyplot:

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

Ответ 5

Двумерная гауссовская матрица ядра может быть вычислена с помощью широковещательной передачи numpy,

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

Для небольших размеров ядра это должно быть достаточно быстро.

Примечание: это упрощает изменение параметра сигмы относительно принятого ответа.

Ответ 6

linalg.norm принимает параметр axis. При небольшом эксперименте я обнаружил, что могу рассчитать норму для всех комбинаций строк с

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

Он расширяет x в 3d-массив всех различий и принимает норму в последнем измерении.

Поэтому я могу применить это к вашему коду, добавив параметр axis к вашему 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)

производит

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

Ответ 7

На основе ответа Тедди Хартанто. Вы можете просто вычислить свои собственные одномерные гауссовские функции, а затем использовать np.outer для вычисления двухмерного. Очень быстрый и эффективный способ.

С помощью приведенного ниже кода вы также можете использовать разные Sigmas для каждого измерения

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

Ответ 8

Я пытался использовать только NumPy. Вот код

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=np.sqrt((i-center)**2+(j-center)**2)
          kernel[i,j]=np.exp(-(diff**2)/(2*sigma**2))
    return kernel/np.sum(kernel)

Вы можете визуализировать результат, используя:

plt.imshow(get_gauss_kernel(5,1))