numpy中二维数组上的矢量化移动窗口

Question

我正在对二维数组中恒定大小的移动窗口应用操作。是否有一个有效的类似矢量化的操作，我可以在不循环 Python 的情况下实现这一点？我目前的结构看起来像这样

 for i in range(1,xmax-1):
     for j in range(1,ymax-1):
        out[i][j] = f(in[i][j],in[i+1][j],in[i-1][j],in[i][j+1],in[i][j-1],...)

eat 在这个问题中留下的 cmets 暗示了将这个操作向量化的可能性，但没有进一步的细节 vectorized indexing/slicing in numpy/scipy?

Answer 1

您可以按照here、here 和here 的说明使用滚动窗口技术，但适用于二维数组。

NumPy中二维滚动窗口的源码：

# Rolling window for 2D arrays in NumPy
import numpy as np

def rolling_window(a, shape):  # rolling window for 2D array
    s = (a.shape[0] - shape[0] + 1,) + (a.shape[1] - shape[1] + 1,) + shape
    strides = a.strides + a.strides
    return np.lib.stride_tricks.as_strided(a, shape=s, strides=strides)

a = np.array([[0,  1,  2,  3,  4,  5],
              [6,  7,  8,  9, 10,  11],
              [12, 13, 14, 15, 7,   8],
              [18, 19, 20, 21, 13, 14],
              [24, 25, 26, 27, 19, 20],
              [30, 31, 32, 33, 34, 35]], dtype=np.int)
b = np.arange(36, dtype=np.float).reshape(6,6)
present = np.array([[7,8],[13,14],[19,20]], dtype=np.int)
absent  = np.array([[7,8],[42,14],[19,20]], dtype=np.int)

found = np.all(np.all(rolling_window(a, present.shape) == present, axis=2), axis=2)
print(np.transpose(found.nonzero()))
found = np.all(np.all(rolling_window(b, present.shape) == present, axis=2), axis=2)
print(np.transpose(found.nonzero()))
found = np.all(np.all(rolling_window(a, absent.shape) == absent, axis=2), axis=2)
print(np.transpose(found.nonzero()))

数组present在数组a中出现了两次，分别在[1,1]和[2,4]上。

我的 CoLab 笔记本"Rolling window on NumPy arrays without for loops" 中有更多示例。

Answer 2

如果你能表达功能

f(in[i][j],in[i+1][j],in[i-1][j],in[i][j+1],in[i][j-1],…)

作为线性运算符，您可以使用 scipy 的 signal.convolve2d 函数来做到这一点。例如，假设您有一个 50x50 数组 A，并且您想计算第二个数组 B，其中每个元素 b[ij] 是数组 A 中 a[i,j], a[(i-1),j], a[i,(j-1)], a[(i-1),(j-1)] 的平均值。您可以这样做：

A = # your first array
B = numpy.ones((2,2))/4
C = scipy.signal.convolve2d(A,B, 'valid')

当执行卷积时，数组 B “滑动”穿过 A，将相应的元素相乘，然后将结果相加。由于边框效果，在使用生成的数组 C 时必须小心。这里，C 的形状为 49x49，因为 convolve2d 中的 'valid' 参数，丢弃第一行和第一列，因为它们包含边框效果。如果你想要一个 50x50 的数组而不丢弃，你可以将该参数换成 'same'

编辑：如果您能告诉我更多关于您需要的函数的信息，我可以更具体地帮助您将其转换为用于进行 2D 卷积的数组。

希望有帮助！