高效密度函数计算

Question

我有一个大的numpy数组形式的图像（opencv将其作为3个uint8值的2D数组返回），想要计算每个像素的高斯核密度和，即（SO中仍没有LaTeX支持吗？）： 其中N个不同的核具有指定权重w、均值和对角协方差矩阵。

所以基本上我想要一个函数compute_densities(image, kernels) -> numpy array of floats。在Python中，最有效的方法是什么？如果scipy中没有库函数进行此操作，我会感到惊讶，但我很久以前就学过统计学，所以我有点困惑于文档的细节。

基本上我想要以下内容，只是比naive 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文件，它会返回一个包含三个颜色通道的三元素uint8结构的二维数组（高x宽）。

Kernels是一个namedtuple('Kernel', 'mean covariance weight')，其中mean是一个向量，covariance是一个除对角线外其余元素均为零的3x3矩阵，而weight是一个浮点数，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)

- Voo

---

请查看scipy.ndimage（http://docs.scipy.org/doc/scipy/reference/ndimage.html）。构建您自己的卷积核，然后使用convolve方法将它们与数组进行卷积。 - dnf0

添加悬赏，希望有人能提供一个示例，说明我应该如何使用卷积和 co 来实现这个函数。 - Voo

Voo，请问你能否在示例中添加一些样本值和对compute_probabilities()进行一个示例调用吗？（以及它会产生什么结果）我会尝试完成这个任务，但输入类型并不完全明显。我认为协方差是一个KxK数组，其中K是img.shape[2]，是正确的吗？ - Alex I

@Alex，我们来了，希望这样更清晰明了。 - Voo

@Voo，请看下面的答案，我认为它可行。 - Alex I

Answer 1

编辑

我验证了这段代码生成的结果与原始代码相同：

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,1）变为（3,）。

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 N x K 的数组和一个 M x K x N 的数组进行点积将产生结果 M x N x M x N，在一维和多维数据上，点积的行为是不同的。因此，我们只需进行元素逐个相乘，然后沿最后一个维度求和即可。

实际上，问题中的“高斯核总和”部分让我感到困惑。所请求的计算方式是，对于每个输出像素，该值仅依赖于同一像素的输入值，而不依赖于相邻像素的值。因此，这与高斯模糊毫不相似（它将使用卷积），而只是对每个像素单独执行的计算。

P.S. 1 / covariance 会有问题。您确定不要使用 np.linalg.inv(covariance) 吗？

旧答案

听起来你想要其中之一：

scipy.signal.convolve2d

scipy.ndimage.filters.convolve

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

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

- Alex I

---

一张图像有几个高斯（大约十几个，在“编译时”已知）。我不确定我的内核是否可分离，但为了澄清，我已经添加了朴素的Python代码。 - Voo

玩了一下卷积，但不知道如何使用它来实现该函数。 - Voo

Answer 2

据我所理解，您的目标是针对多元正态分布混合模型评估每个图像像素值的概率密度。

[当前（2013-11-20）您问题中的代码以及@Alex I的回答存在一个错误 - 上述方程式中的“||”实际上表示行列式而不是向量范数 - 例如请参见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）。截至0.13版本，SciPy中最接近开箱即用的功能是scipy.stats中的核密度估计实现（请参见here），它定义了一个非常类似的分布，其中每个内核的协方差矩阵都是相同的-因此不适合您的问题。

- rroowwllaanndd

---

哦，我从来没有见过在行列式中使用“||”，而且对于我的小例子选择，结果与正确解决方案非常相似，以至于我甚至没有注意到，这本身就应该赢得您的奖励;) 太棒了 - 我必须更多地使用numpy，但是重新塑形和步幅似乎是解决许多非直观情况的最佳选择，既高效又容易。 - Voo

Answer 3

使用Python获得更好的性能的方法是不使用Python。有许多包可以与Python语法一起使用，但是使用C或C ++后端。NumPy本身就是这样做的。您的问题似乎非常适合使用Cython或numexpr之类的工具。这两个链接都向您展示了如何在NumPy向量上使用这两种系统的内核。 编辑：我希望我的其中一个投票者能让我知道我错了。如果找不到预制函数，我建议一种方法。如果您知道比Cython或numexpr更高效的方法（即以Python语法编写C的方法），那么我很想听听它。

- Adam

---

我认为已经有一些统计软件可以做到这一点，因为这似乎并不是什么罕见的事情。但是如果必须的话，我会用Cython来编写它，只是感觉应该已经有类似的东西存在！ - Voo

也许你的开头语句就是被投票者所阅读的全部内容。 - Geoff

@Geoff 你可能是对的。我可能需要找到一种更为委婉的说法。尽管那个陈述是令人遗憾的真实。 - Adam