Python 如何在numpy中有效地计算高斯核矩阵？
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/29731726/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
How to calculate a Gaussian kernel matrix efficiently in numpy?
提问by hidemyname
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.
这是我目前的方式。有什么办法可以使用矩阵运算来做到这一点吗？X 是数据点。
采纳答案by FuzzyDuck
Do you want to use the Gaussian kernel for e.g. image smoothing? If so, there's a function gaussian_filter()in scipy:
你想使用高斯核来进行图像平滑吗？如果是这样，gaussian_filter()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.htmfor an example.
这应该有效 - 虽然它仍然不是 100% 准确，但它试图解释网格每个单元格内的概率质量。我认为在每个单元格的中点使用概率密度稍微不太准确，尤其是对于小内核。有关示例，请参见https://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm。
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:
通过链接在图 3 中的示例上对其进行测试：
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 wrongThe 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
回答by hpaulj
linalg.normtakes an axisparameter. With a little experimentation I found I could calculate the norm for all combinations of rows with
linalg.norm接受一个axis参数。通过一些实验，我发现我可以计算所有行组合的范数
np.linalg.norm(x[None,:,:]-x[:,None,:],axis=2)
It expands xinto a 3d array of all differences, and takes the norm on the last dimension.
它扩展x为所有差异的 3d 数组，并在最后一个维度上取范数。
So I can apply this to your code by adding the axisparameter to your Gaussian:
因此，我可以通过将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)
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]])
回答by Nils Werner
You may simply gaussian-filter a simple 2D dirac function, the result is then the filter function that was being used:
您可以简单地高斯过滤一个简单的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)
回答by clemisch
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)
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)
Edit: Changed arange to linspace to handle even side lengths
编辑：将 arange 更改为 linspace 以处理均匀的边长
回答by rth
A 2D gaussian kernel matrix can be computed with numpy broadcasting,
可以使用 numpy 广播计算 2D 高斯核矩阵，
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.
注意：这使得相对于接受的答案更容易更改 sigma 参数。
回答by Teddy Hartanto
I'm trying to improve on FuzzyDuck's answerhere. I think this approach is shorter and easier to understand. Here I'm using signal.scipy.gaussianto get the 2D gaussian kernel.
我正在尝试改进FuzzyDuck 的回答。我认为这种方法更短也更容易理解。在这里，我使用signal.scipy.gaussian2D 高斯内核。
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:
绘制它使用matplotlib.pyplot：
import matplotlib.pyplot as plt
plt.imshow(gkern(21), interpolation='none')
回答by Dodo
Building up on Teddy Hartanto's answer. You can just calculate your own one dimensional Gaussian functions and then use np.outerto calculate the two dimensional one. Very fast and efficient way.
以 Teddy Hartanto 的回答为基础。您可以计算自己的一维高斯函数，然后使用它np.outer来计算二维一维。非常快速有效的方式。
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
回答by Suraj Singh
I tried using numpy only. Here is the code
我尝试只使用 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)
You can visualise the result using:
您可以使用以下方法可视化结果：
plt.imshow(get_gauss_kernel(5,1))
回答by Vinoj John Hosan
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)
最大在 : (2, 2)
kernel shape (10, 10)
内核形状 (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)
最大 : (40, 25)
kernel shape (50, 100)
内核形状 (50, 100)
回答by Prasad Raghavendra
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.
例如，这是 3x5 补丁的单行函数。
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.
您可以在此处阅读有关 scipy 高斯的更多信息。
