如何在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)) 

这是我目前的做法。有什么办法可以使用矩阵操作来做到这一点？ X是数据点。

Answer 1

您是否想使用高斯内核来处理图像平滑？如果是这样，有一个功能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 

输入：

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

输出：

Answer 2

linalg.norm需要一个axis参数。随着一点点的实验，我发现我可以计算出中行的所有组合的规范与

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

它扩大到x所有差异的三维阵列，并采取在最后一维的常态。

所以我可以通过添加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]]) 

匹配：

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

Answer 3

你可能只是高斯过滤一个简单的2D dirac function，结果就是过滤器f这是正在使用的油膏：

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 4

我自己用公认的答案对我的图像处理，但我觉得它（以及其他的答案）过于依赖于其他模块。此外，接受的答案有时会产生内核，并且最终会有很多零条目。

因此，这里是我的紧凑型解决方案：

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) 

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

对于小尺寸的内核，这应该是相当快的。

注：这使得相对于接受的答案更容易改变西格玛参数。

Answer 6

我想在这里改进FuzzyDuck's answer。我认为这种方法更短，更易于理解。在这里，我使用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') 

Answer 7

上泰迪Hartanto的答案构建。您可以计算您自己的一维高斯函数，然后使用np.outer来计算二维函数。非常快速和有效的方式。

有了下面的代码，您还可以使用不同的Sigma的每一个维度

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