高效的密度函数计算

Question

我有一个numpy数组形式的大图像(opencv将其作为3个uint8值的2d数组返回)并且想要为每个像素计算高斯内核的总和,即(在那里仍然没有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结构的2d数组(高x宽).

内核是a namedtuple('Kernel', 'mean covariance weight'),mean是一个向量,协方差是一个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

编辑

我确认这会产生与原始代码相同的结果:

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,点在一维和多维数据上表现不同.所以我们只做一个逐元素乘法,然后沿最后一个维度求和.

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

PS 1 / covariance是有问题的.你确定你不想要np.linalg.inv(covariance)吗？

老答复

这听起来像你想要的是其中之一:

scipy.signal.convolve2d

scipy.ndimage.filters.convolve

问题有点令人困惑,你试图计算一堆用不同的高斯卷积的图像,或者用一个高斯的数量卷积的单个图像？你的内核是可分的吗？(如果是,使用两个卷积Mx1和1xN而不是一个MxN)你使用的scipy函数在任何情况下都是相同的.

当然也你想要用的组合预先计算的内核numpy.random.normal和meshgrid.