python - 高效的密度函数计算

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 个颜色 channel 的 3 uint8 结构的二维数组(高 x 宽)。

内核是一个 namedtuple('Kernel', 'mean covariance weight') ，均值是一个向量，协方差是一个3x3除对角线以外的所有内容均为零且权重为 float 的矩阵 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

编辑

我验证了这会产生与原始代码相同的结果:

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 N x K 点缀数组 M x K x N 会产生结果 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 的组合来预先计算您的内核。