高效的密度函数计算

Question

我有一个 numpy 数组形式的大图像（opencv 将它作为 3 个 uint8 值的二维数组返回）并且想要计算每个像素的高斯内核的总和，即（SO 中仍然没有 LaTeX 支持吗？ ):

对于具有指定权重 w、均值和对角协方差矩阵的 N 个不同内核。

所以基本上我想要一个函数compute_densities(image, kernels) -> numpy array of floats。在 python 中有效地做到这一点的最佳方法是什么？如果 scipy 中还没有为此的库函数，我会感到惊讶，但我很久以前在 uni 有统计数据，所以我对文档的细节有点困惑..

基本上我想要以下，只是比天真的 python (2pi^{-3/2} 被忽略，因为它是一个常数因素，对我来说并不重要，因为我只对概率之间的比率感兴趣)

def compute_probabilities(img, kernels):
    np.seterr(divide='ignore') # 1 / covariance logs an error otherwise
    result = np.zeros((img.shape[0], img.shape[1]))
    for row_pos, row_val in enumerate(img):
        for col_pos, val in enumerate(row_val):
            prob = 0.0
            for kernel in kernels:
                mean, covariance, weight = kernel
                val_sub_mu = np.array([val]).T - mean
                cov_inv = np.where(covariance != 0, 1 / covariance, 0)
                tmp = val_sub_mu.T.dot(cov_inv).dot(val_sub_mu)
                prob += weight / np.sqrt(np.linalg.norm(covariance)) * \
                        math.exp(-0.5 * tmp)
            result[row_pos][col_pos] = prob
    np.seterr(divide='warn')
    return result

输入：cv2.imread 在一些 jpg 上，它给出了一个包含 3 个颜色通道的 3 uint8 结构的二维数组（高 x 宽）。

内核是一个namedtuple('Kernel', 'mean covariance weight')，平均值是一个向量，协方差是一个3x3 矩阵，除了对角线为零，权重是一个浮点数0 < weight < 1。为简单起见，我只指定对角线，然后将其转换为 3x3 矩阵：（表示不是一成不变的，我不在乎它是如何表示的，所以可以随意更改所有这些）：

some_kernels = [
   Kernel(np.array([(73.53, 29.94, 17.76)]), np.array([(765.40, 121.44, 112.80)]), 0.0294),
   ...
]

def fixup_kernels(kernels):
    new_kernels = []
    for kernel in kernels:
        cov = np.zeros((3, 3))
        for pos, c in enumerate(kernel.covariance[0]):
            cov[pos][pos] = c
        new_kernels.append(Kernel(kernel.mean.T, cov, kernel.weight))
    return new_kernels

 some_kernels = fixup_kernels(some_kernels)
 img = cv2.imread("something.jpg")
 result = compute_probabalities(img, some_kernels)

Answer 1

据我了解，您的目标是评估每个图像像素值相对于多元正态分布混合的概率密度。

[目前 (2013-11-20) 您的问题和@Alex I 的回答中的代码有错误 - 上述等式中围绕\Sigma 的| | 实际上表示行列式 而不是向量范数 - 参见例如here。在对角协方差的情况下，行列式只是对角元素的乘积。]

可以非常高效地实现密度计算 就numpy数组操作而言。以下实现利用 问题中协方差矩阵的球形（即对角线）性质：

def compute_probabilities_faster(img, kernels):
  means, covs, weights = map(np.dstack, zip(*kernels)) 
  pixels_as_rows = img.reshape((-1, 3, 1))
  responses = np.exp(-0.5 * ((pixels_as_rows - means) ** 2 / covs).sum(axis=1))
  factors = 1. / np.sqrt(covs.prod(axis=1) * ((2 * np.pi) ** 3))
  return np.sum(responses * factors * weights, axis=2).reshape(img.shape[:2])

这个函数直接在内核上运行，就像它们最初一样 代表，即。 没有由您的fixup_kernels 函数修改。当规范化因子(2 * np.pi) ** 3 被移除（并且对linalg.norm 的调用被替换为linalg.det）时，此函数与您的代码输出相匹配（足以满足np.allclose）。

SciPy 中最接近的开箱即用功能（截至 0.13）是 scipy.stats 中内核密度估计的实现（参见here），它定义了一个非常相似的 每个内核的协方差矩阵相同的分布 - 因此不适合您的问题。

Answer 2

编辑

我证实这会产生与原始代码相同的结果：

def compute_probabilities_fast(img, kernels):
    np.seterr(divide='ignore')
    result = np.zeros((img.shape[0], img.shape[1]))
    for kernel in kernels:
        mean, covariance, weight = kernel
        cov_inv = np.where(covariance != 0, 1 / covariance, 0)
        mean = mean[:,0]
        img_sub_mu = img - mean
        img_tmp = np.sum( img_sub_mu.dot(cov_inv) * img_sub_mu, axis=2 )
        result += (weight / np.sqrt(np.linalg.norm(covariance))) * np.exp(-0.5 * img_tmp)
    return result

解释：

mean[:,0] 将形状简化为 (3,) 而不是 (3,1)。

img - mean 广播到整个图像并从每个像素中减去均值。

img_sub_mu.dot(cov_inv) 大致相当于val_sub_mu.T.dot(cov_inv)。

np.sum( ... * img_sub_mu, axis=2 ) 大致相当于.dot(val_sub_mu)。但是，不能使用点，因为这样做会增加额外的尺寸。例如，用数组 M x K x N 点缀的数组 M x N x K 将产生结果 M x N x M x N，点在一维和多维数据上的行为不同。所以我们只做一个元素乘法，然后沿着最后一个维度求和。

实际上，问题中的“高斯核之和”让我感到困惑。所要求的是计算，对于每个输出像素，其值仅取决于同一像素的输入值，而不取决于相邻像素的值。因此，这与高斯模糊（将使用卷积）完全不同，它只是对每个像素单独执行的计算。

附： 1 / covariance 有问题。你确定不想要np.linalg.inv(covariance) 吗？

老答案

听起来你想要的是以下之一：

scipy.signal.convolve2d

scipy.ndimage.filters.convolve

这个问题有点令人困惑，你是想计算一堆用不同高斯卷积的图像，还是用高斯总和卷积的单个图像？你的内核是可分离的吗？ （如果是，请使用两个卷积 Mx1 和 1xN 而不是一个 MxN）您使用的 scipy 函数在任何情况下都是相同的。

当然，您还想使用numpy.random.normal 和meshgrid 的组合预先计算您的内核。

Answer 3

从 Python 获得性能的方法是不使用 Python。

那里有许多使用 Python 语法的包，但随后使用 C 或 C++ 后端。 NumPy 本身就是这样做的。您的问题似乎是为Cython 或numexpr 等量身定制的。这两个链接都向您展示了如何将任一系统用于 NumPy 向量上的内核。

编辑：我希望我的一位反对者告诉我我错了。如果找不到预制功能，我建议采取一种方法。如果您知道一种比 Cython 或 numexpr（即用 Python 语法编写 C 的方法）具有更高性能的方法，那么我很想听听。