更快地计算图像(M,N)和模板(3,3)之间的平方差之和以进行模板匹配？

Question

我正在实现这里概述的纹理合成算法.为此,我需要计算平方差的和,这是一个度量,用于估计template不同位置之间的误差image.我的工作实施缓慢如下:

total_weight = valid_mask.sum()  
for i in xrange(input_image.shape[0]):  
    for j in xrange(input_image.shape[1]):  
        sample = image[i:i + window, j:j + window]  
        dist = (template - sample) ** 2  
        ssd[i, j] = (dist * valid_mask).sum() / total_weight  

这里,total_weight只是为了规范化.有些像素的强度未知,因此我用它valid_mask来掩盖它们.这个嵌套循环位于2个循环内部,因此这是4个嵌套循环,这显然是性能杀手!

有没有办法让我在NumPy或Python中更快,这个嵌套循环的替代品？矢量化是可能的吗？我需要在上班(3, 3)的部分image用(3,3)的template.

我随后将在Cython中实现这一点,所以我可以更快地使用NumPy工作,更好的是.

你可以在这里找到完整的代码.第62-67行引用了这里.

谢谢,
Chintak

Answer 1

这基本上是对Warren Weckesser答案的改进.前进的方法显然是使用原始数组的多维窗口视图,但是您希望保持该视图不会触发副本.如果你扩展你的sum((a-b)**2),你可以把它变成sum(a**2) + sum(b**2) - 2*sum(a*b),然后你可以用线性代数运算符执行这个乘法 - 然后 - 减少 - 和 - 运算,在性能和内存使用方面有显着改进:

def sumsqdiff3(input_image, template):
    window_size = template.shape
    y = as_strided(input_image,
                    shape=(input_image.shape[0] - window_size[0] + 1,
                           input_image.shape[1] - window_size[1] + 1,) +
                          window_size,
                    strides=input_image.strides * 2)
    ssd = np.einsum('ijkl,kl->ij', y, template)
    ssd *= - 2
    ssd += np.einsum('ijkl, ijkl->ij', y, y)
    ssd += np.einsum('ij, ij', template, template)

    return ssd

In [288]: img = np.random.rand(500, 500)

In [289]: template = np.random.rand(3, 3)

In [290]: %timeit a = sumsqdiff2(img, template) # Warren's function
10 loops, best of 3: 59.4 ms per loop

In [291]: %timeit b = sumsqdiff3(img, template)
100 loops, best of 3: 18.2 ms per loop

In [292]: np.allclose(a, b)
Out[292]: True

我valid_mask故意将参数排除在外,因为我不完全理解你将如何使用它.原则上,只需将相应的值归零template和/或input_image应该执行相同的技巧.

Answer 2

你可以as_strided结合numpy的广播功能做一些惊人的事情.以下是您的两个版本的功能:

import numpy as np
from numpy.lib.stride_tricks import as_strided


def sumsqdiff(input_image, template, valid_mask=None):
    if valid_mask is None:
        valid_mask = np.ones_like(template)
    total_weight = valid_mask.sum()
    window_size = template.shape
    ssd = np.empty((input_image.shape[0] - window_size[0] + 1,
                    input_image.shape[1] - window_size[1] + 1))
    for i in xrange(ssd.shape[0]):  
        for j in xrange(ssd.shape[1]):  
            sample = input_image[i:i + window_size[0], j:j + window_size[1]]  
            dist = (template - sample) ** 2  
            ssd[i, j] = (dist * valid_mask).sum()
    return ssd


def sumsqdiff2(input_image, template, valid_mask=None):
    if valid_mask is None:
        valid_mask = np.ones_like(template)
    total_weight = valid_mask.sum()
    window_size = template.shape

    # Create a 4-D array y, such that y[i,j,:,:] is the 2-D window
    #     input_image[i:i+window_size[0], j:j+window_size[1]]
    y = as_strided(input_image,
                    shape=(input_image.shape[0] - window_size[0] + 1,
                           input_image.shape[1] - window_size[1] + 1,) +
                          window_size,
                    strides=input_image.strides * 2)

    # Compute the sum of squared differences using broadcasting.
    ssd = ((y - template) ** 2 * valid_mask).sum(axis=-1).sum(axis=-1)
    return ssd

这是一个比较它们的ipython会话.

我将用于演示的模板:

In [72]: template
Out[72]: 
array([[-1,  1, -1],
       [ 1,  2,  1],
       [-1,  1, -1]])

一个小输入,所以我们可以检查结果:

In [73]: x
Out[73]: 
array([[  0.,   1.,   2.,   3.,   4.,   5.,   6.],
       [  7.,   8.,   9.,  10.,  11.,  12.,  13.],
       [ 14.,  15.,  16.,  17.,  18.,  19.,  20.],
       [ 21.,  22.,  23.,  24.,  25.,  26.,  27.],
       [ 28.,  29.,  30.,  31.,  32.,  33.,  34.]])

应用这两个函数x并检查我们得到相同的结果:

In [74]: sumsqdiff(x, template)
Out[74]: 
array([[  856.,  1005.,  1172.,  1357.,  1560.],
       [ 2277.,  2552.,  2845.,  3156.,  3485.],
       [ 4580.,  4981.,  5400.,  5837.,  6292.]])

In [75]: sumsqdiff2(x, template)
Out[75]: 
array([[  856.,  1005.,  1172.,  1357.,  1560.],
       [ 2277.,  2552.,  2845.,  3156.,  3485.],
       [ 4580.,  4981.,  5400.,  5837.,  6292.]])

现在做一个更大的输入"图像":

In [76]: z = np.random.randn(500, 500)

并检查性能:

In [77]: %timeit sumsqdiff(z, template)
1 loops, best of 3: 3.55 s per loop

In [78]: %timeit sumsqdiff2(z, template)
10 loops, best of 3: 33 ms per loop

不是太寒酸.:)

两个缺点:

• 计算sumsqdiff2将生成一个临时数组,对于3x3模板,将是9倍的大小input_image.(一般来说,它template.size的大小将是input_image.)
• 在您对代码进行Cython化时,这些"跨步技巧"对您无济于事.转换为Cython时,你经常会在使用numpy进行矢量化时放回你摆脱的循环中.