如何在Python中高效计算两个高斯分布的热力图？

Question

我正在尝试生成一个热力图，其中像素值由两个独立的二维高斯分布控制。让它们分别为Kernel1（muX1，muY1，sigmaX1，sigmaY1）和Kernel2（muX2，muY2，sigmaX2，sigmaY2）。更具体地说，每个内核的长度是其标准偏差的三倍。第一个内核具有sigmaX1 = sigmaY1，而第二个内核具有sigmaX2 < sigmaY2。两个内核的协方差矩阵都是对角线矩阵（X和Y是独立的）。Kernel1通常完全包含在Kernel2中。

我尝试了以下两种方法，但结果都不理想。能否给我一些建议？

方法1：

遍历地图上的所有像素值对（i，j），并计算给定I（i，j）= P（i，j | Kernel1，Kernel2）的I（i，j）的值= 1-（1-P（i，j | Kernel1））*（1-P（i，j | Kernel2））。然后我得到了以下结果，在平滑性方面效果良好。但它在我的电脑上运行需要10秒钟，速度太慢了。

代码：

def genDensityBox(self, height, width, muY1, muX1, muY2, muX2, sigmaK1, sigmaY2, sigmaX2):
    densityBox = np.zeros((height, width))
    for y in range(height):
        for x in range(width):
            densityBox[y, x] += 1. - (1. - multivariateNormal(y, x, muY1, muX1, sigmaK1, sigmaK1)) * (1. - multivariateNormal(y, x, muY2, muX2, sigmaY2, sigmaX2))
    return densityBox

def multivariateNormal(y, x, muY, muX, sigmaY, sigmaX):
    return norm.pdf(y, loc=muY, scale=sigmaY) * norm.pdf(x, loc=muX, scale=sigmaX)

第二种方法：

分别生成两个与两个核对应的图像，然后使用特定的alpha值将它们混合在一起。每个图像都是通过取两个一维高斯滤波器的外积来生成的。然后我得到了以下非常粗糙的结果。但这种方法的优点是由于在两个向量之间使用外积，所以速度非常快。

由于第一种方法很慢而第二种方法很粗糙，我正在尝试找到一种新的方法，既可以实现良好的平滑度，又可以同时具有低时间复杂度。有人能给我一些帮助吗？

谢谢！

对于第二种方法，可以像这里提到的那样轻松生成2D高斯地图：here：

def gkern(self, sigmaY, sigmaX, yKernelLen, xKernelLen, nsigma=3):
    """Returns a 2D Gaussian kernel array."""
    yInterval = (2*nsigma+1.)/(yKernelLen)
    yRow = np.linspace(-nsigma-yInterval/2.,nsigma+yInterval/2.,yKernelLen + 1)
    kernelY = np.diff(st.norm.cdf(yRow, 0, sigmaY))
    xInterval = (2*nsigma+1.)/(xKernelLen)
    xRow = np.linspace(-nsigma-xInterval/2.,nsigma+xInterval/2.,xKernelLen + 1)
    kernelX = np.diff(st.norm.cdf(xRow, 0, sigmaX))    
    kernelRaw = np.sqrt(np.outer(kernelY, kernelX))
    kernel = kernelRaw / (kernelRaw.sum())
    return kernel

- Louis Kuang

Answer 1

你的方法很好，除了不应该循环使用 norm.pdf，而应该将所有需要计算核函数的值都推送到一起，然后将输出重新整形为图像的所需形状。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# create 2 kernels
m1 = (-1,-1)
s1 = np.eye(2)
k1 = multivariate_normal(mean=m1, cov=s1)

m2 = (1,1)
s2 = np.eye(2)
k2 = multivariate_normal(mean=m2, cov=s2)

# create a grid of (x,y) coordinates at which to evaluate the kernels
xlim = (-3, 3)
ylim = (-3, 3)
xres = 100
yres = 100

x = np.linspace(xlim[0], xlim[1], xres)
y = np.linspace(ylim[0], ylim[1], yres)
xx, yy = np.meshgrid(x,y)

# evaluate kernels at grid points
xxyy = np.c_[xx.ravel(), yy.ravel()]
zz = k1.pdf(xxyy) + k2.pdf(xxyy)

# reshape and plot image
img = zz.reshape((xres,yres))
plt.imshow(img); plt.show()

这种方法不应该花费太长时间：

In [26]: %timeit zz = k1.pdf(xxyy) + k2.pdf(xxyy)
1000 loops, best of 3: 1.16 ms per loop

- Paul Brodersen

Answer 2

根据Paul的答案，我制作了一个函数，可以根据高斯分布的中心生成热力图（这可能对其他人有帮助）：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

def points_to_gaussian_heatmap(centers, height, width, scale):
    gaussians = []
    for y,x in centers:
        s = np.eye(2)*scale
        g = multivariate_normal(mean=(x,y), cov=s)
        gaussians.append(g)

    # create a grid of (x,y) coordinates at which to evaluate the kernels
    x = np.arange(0, width)
    y = np.arange(0, height)
    xx, yy = np.meshgrid(x,y)
    xxyy = np.stack([xx.ravel(), yy.ravel()]).T
    
    # evaluate kernels at grid points
    zz = sum(g.pdf(xxyy) for g in gaussians)

    img = zz.reshape((height,width))
    return img

W = 800  # width of heatmap
H = 400  # height of heatmap
SCALE = 64  # increase scale to make larger gaussians
CENTERS = [(100,100), 
           (100,300), 
           (300,100)] # center points of the gaussians

img = points_to_gaussian_heatmap(CENTERS, H, W, SCALE)

plt.imshow(img); plt.show()

- Ismael EL ATIFI

---

我认为“Y”轴被倒置索引了。 - undefined

这只是因为pyplot.imshow在默认情况下处理原点的方式。可以通过指定origin参数为plt.imshow(img, origin='lower')来修复这个问题。 - undefined