7

我正在实现此处概述的纹理合成算法。为此，我需要计算平方差之和，这是一个用于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.template

我随后将在 Cython 中实现它，所以我可以使用 NumPy 越快地让它工作，它就越好。

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

谢谢，
钦塔克

4

4 回答 4

10

这基本上是对 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应该做同样的事情。

于 2013-07-26T16:01:55.073 回答
6

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 模板，该数组的大小将是input_image. （通常它template.size的大小是 的倍数input_image。）
  • 当您对代码进行 Cythonize 时，这些“跨步技巧”对您没有帮助。转换为 Cython 时，您通常最终会回到使用 numpy 进行矢量化时摆脱的循环。
于 2013-07-26T13:51:20.017 回答
1

如果您重新排列算法以逐行执行计算，那么检查它的执行情况可能是值得的。这个想法是，如果连续读取内存，您可能会更好地使用 CPU 缓存。

伪代码：

for template_row in template:
  for row in image:
    for col in image:
      # find distance template_row to sample_row
      # add sum to ssd[row - template_row, col]

实际代码（在沃伦之后）：

def sumsqdiffr(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.zeros((input_image.shape[0] - window_size[0] + 1,
                    input_image.shape[1] - window_size[1] + 1))

    for tr in xrange(template.shape[0]):
        for i in xrange(tr, ssd.shape[0] + tr):
            for j in xrange(ssd.shape[1]):  
                sample = input_image[i, j:j + window_size[1]]  
                dist = (template[tr] - sample) ** 2  
                ssd[i - tr, j] += (dist * valid_mask[tr]).sum()
    return ssd

它比原始实现慢两倍多。

（如果有人想告诉我整个想法是错误的还是导致这种情况的原因，我很乐意从中获得一些理解）

于 2013-07-26T13:31:02.923 回答
0

我认为您在实现算法方面做得很好。矢量化是一种选择，但我鼓励您使用Numba优化编译器，它将 Python 语法编译为机器代码。Numba 效果令人印象深刻。Numba vs. Cython：Take 2是对 Numba 的非常简短的介绍和性能比较。

于 2013-07-26T12:59:28.163 回答