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

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

Это мой нынешний путь. Есть ли способ использовать матричную операцию для этого? X - это точки данных.

7 ответов

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

В качестве альтернативы это должно работать:

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

Input:

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

Вывод:   


Вы можете просто 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)


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

Итак, вот мое компактное решение:

import numpy as np
def gkern(l=5, sig=1.):
 """
 creates gaussian kernel with side length l and a sigma of sig
 """
 ax = np.arange(-l // 2 + 1., l // 2 + 1.)
 xx, yy = np.meshgrid(ax, ax)
 kernel = np.exp(-(xx**2 + yy**2) / (2. * sig**2))
 return kernel / np.sum(kernel)


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


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

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:

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


Двумерная гауссовская матрица ядра может быть вычислена с помощью широковещательной передачи 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()

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

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


На основе ответа Тедди Хартанто. Вы можете просто вычислить свои собственные одномерные гауссовские функции, а затем использовать 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

licensed under cc by-sa 3.0 with attribution.