Python - 矢量化滑动窗口

Question

我正在尝试矢量化滑动窗口操作。对于一维情况，一个有用的例子可以是：

x= vstack((np.array([range(10)]),np.array([range(10)])))

x[1,:]=np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]+1],x[1,:])

索引

x[1,:]=np.where((x[0,:]<2)&(x[0,:]>0),x[1,x[0,:]+1],x[1,:])
IndexError: index (10) out of range (0<=index<9) in dimension 1

奇怪的是，对于 n-1 值，我不会得到这个错误，这意味着索引小于 0。它似乎并不介意：

x[1,:]=np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-1],x[1,:])

print(x)

[[0 1 2 3 4 5 6 7 8 9]
 [0 0 1 2 3 5 6 7 8 9]]

这有什么问题吗？我的方法完全错误吗？任何 cmets 将不胜感激。

编辑：

这就是我想要实现的，我将一个矩阵展平为一个 numpy 数组，我想在该数组上计算每个单元格的 6x6 邻域的平均值：

matriz = np.array([[1,2,3,4,5],
   [6,5,4,3,2],
   [1,1,2,2,3],
   [3,3,2,2,1],
   [3,2,1,3,2],
   [1,2,3,1,2]])

# matrix to vector
vector2 = ndarray.flatten(matriz)

ncols = int(shape(matriz)[1])
nrows = int(shape(matriz)[0])

vector = np.zeros(nrows*ncols,dtype='float64')


# Interior pixels
if ( (i % ncols) != 0 and (i+1) % ncols != 0 and i>ncols and i<ncols*(nrows-1)):

    vector[i] = np.mean(np.array([vector2[i-ncols-1],vector2[i-ncols],vector2[i-ncols+1],vector2[i-1],vector2[i+1],vector2[i+ncols-1],vector2[i+ncols],vector2[i+ncols+1]]))

Answer 1

如果我正确理解了这个问题，您希望在索引周围取所有数字的平均值 1 步，忽略索引。

我已经修补了你的功能，我相信你会这样做：

def original(matriz):

    vector2 = np.ndarray.flatten(matriz)

    nrows, ncols= matriz.shape
    vector = np.zeros(nrows*ncols,dtype='float64')

    # Interior pixels
    for i in range(vector.shape[0]):
        if ( (i % ncols) != 0 and (i+1) % ncols != 0 and i>ncols and i<ncols*(nrows-1)):

            vector[i] = np.mean(np.array([vector2[i-ncols-1],vector2[i-ncols],\
                        vector2[i-ncols+1],vector2[i-1],vector2[i+1],\
                        vector2[i+ncols-1],vector2[i+ncols],vector2[i+ncols+1]]))

我使用切片和视图重写了这个：

def mean_around(arr):
    arr=arr.astype(np.float64)

    out= np.copy(arr[:-2,:-2])  #Top left corner
    out+= arr[:-2,2:]           #Top right corner
    out+= arr[:-2,1:-1]         #Top center
    out+= arr[2:,:-2]           #etc
    out+= arr[2:,2:]
    out+= arr[2:,1:-1]
    out+= arr[1:-1,2:]
    out+= arr[1:-1,:-2]

    out/=8.0    #Divide by # of elements to obtain mean

    cout=np.empty_like(arr)  #Create output array
    cout[1:-1,1:-1]=out      #Fill with out values
    cout[0,:]=0;cout[-1,:]=0;cout[:,0]=0;cout[:,-1]=0 #Set edges equal to zero

    return  cout

使用np.empty_like 然后填充边缘似乎比np.zeros_like 稍微快一些。首先让我们使用您的 matriz 数组仔细检查它们是否提供相同的东西。

print np.allclose(mean_around(matriz),original(matriz))
True

print mean_around(matriz)
[[ 0.     0.     0.     0.     0.   ]
 [ 0.     2.5    2.75   3.125  0.   ]
 [ 0.     3.25   2.75   2.375  0.   ]
 [ 0.     1.875  2.     2.     0.   ]
 [ 0.     2.25   2.25   1.75   0.   ]
 [ 0.     0.     0.     0.     0.   ]]

一些时间安排：

a=np.random.rand(500,500)

print np.allclose(original(a),mean_around(a))
True

%timeit mean_around(a)
100 loops, best of 3: 4.4 ms per loop

%timeit original(a)
1 loops, best of 3: 6.6 s per loop

大约 1500 倍加速。

看起来是个使用 numba 的好地方：

def mean_numba(arr):
    out=np.zeros_like(arr)
    col,rows=arr.shape

    for x in xrange(1,col-1):
        for y in xrange(1,rows-1):
            out[x,y]=(arr[x-1,y+1]+arr[x-1,y]+arr[x-1,y-1]+arr[x,y+1]+\
                      arr[x,y-1]+arr[x+1,y+1]+arr[x+1,y]+arr[x+1,y-1])/8.
    return out

nmean= autojit(mean_numba)

现在让我们与所有提供的方法进行比较。

a=np.random.rand(5000,5000)

%timeit mean_around(a)
1 loops, best of 3: 729 ms per loop

%timeit nmean(a)
10 loops, best of 3: 169 ms per loop

#CT Zhu's answer
%timeit it_mean(a)
1 loops, best of 3: 36.7 s per loop

#Ali_m's answer
%timeit fast_local_mean(a,(3,3))
1 loops, best of 3: 4.7 s per loop

#lmjohns3's answer
%timeit scipy_conv(a)
1 loops, best of 3: 3.72 s per loop

numba up 的 4 倍速度是相当名义上的，这表明 numpy 代码与它所获得的一样好。我提取了其他代码，尽管我确实必须更改@CTZhu 的答案以包含不同的数组大小。

Answer 2

听起来您正在尝试计算 2D 卷积。如果您能够使用scipy，我建议您尝试scipy.signal.convolve2d：

matriz = np.random.randn(10, 10)

# to average a 3x3 neighborhood
kernel = np.ones((3, 3), float)

# to compute the mean, divide by size of neighborhood
kernel /= kernel.sum()

average = scipy.signal.convolve2d(matriz, kernel)

如果将 convolve2d “展开”到其组成循环中，就可以看出计算所有 3x3 邻域的平均值的原因。有效地（并且忽略在源和内核数组边缘发生的事情），它正在计算：

X, Y = kernel.shape
for i in range(matriz.shape[0]):
    for j in range(matriz.shape[1]):
        for ii in range(X):
            for jj in range(Y):
                average[i, j] += kernel[ii, jj] * matriz[i+ii, j+jj]

所以如果你的内核中的每个值都是 1/(1+1+1+1+1+1+1+1+1) == 1/9，你可以将上面的代码重写为：

for i in range(matriz.shape[0]):
    for j in range(matriz.shape[1]):
        average[i, j] = 1./9 * matriz[i:i+X, j:j+Y].sum()

这与在 3x3 区域上计算 matriz 中的值的平均值完全相同，从 i, j 开始。

这样做的一个好处是，您可以通过在内核中适当地设置值来轻松更改与邻域相关的权重。因此，例如，如果你想赋予每个邻域的中心值两倍于其他邻域的权重，你可以这样构建你的内核：

kernel = np.ones((3, 3), float)
kernel[1, 1] = 2.
kernel /= kernel.sum()

卷积码将保持不变，但计算会产生不同类型的平均值（“中心加权”平均值）。这里有很多可能性；希望这为您正在执行的任务提供了一个很好的抽象。

Answer 3

恰好在 Scipy 标准库中有一个函数可以非常快速地计算滑动窗口的平均值。它被称为uniform_filter。您可以使用它来实现邻域均值函数，如下所示：

from scipy.ndimage.filters import uniform_filter
def neighbourhood_average(arr, win=3):
    sums = uniform_filter(arr, win, mode='constant') * (win*win)
    return ((sums - arr) / (win*win - 1))

这将返回一个数组X，其中X[i,j] 是arr 中i,j 的所有邻居的平均值，不包括i,j 本身。请注意，第一列和最后一列以及第一和最后一行受边界条件约束，因此可能对您的应用程序无效（如果需要，您可以使用mode=控制边界规则）。

因为uniform_filter 使用在直接 C 中实现的高效线性时间算法（仅在 arr 的大小上是线性的），它应该很容易胜过任何其他解决方案，尤其是当 win 很大时。

Answer 4

问题出在x[1,x[0,:]+1]，第2轴的索引：x[0,:]+1是[1 2 3 4 5 6 7 8 9 10]，其中索引10大于x的维度。

对于x[1,x[0,:]-1]，第二轴的索引是[-1 0 1 2 3 4 5 6 7 8 9]，你最终会得到[9 0 1 2 3 4 5 6 7 8]，因为9是最后一个元素，索引是-1。倒数第二个元素的索引为-2，以此类推。

对于np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-1],x[1,:]) 和x[0,:]=[0 1 2 3 4 5 6 7 8 9]，本质上是第一个单元格采用x[1,:] 的形式，因为x[0,0] 是0 而x[0,:]<5)&(x[0,:]>0 是False。接下来的四个元素取自x[1,x[0,:]-1]。其余来自x[1,:]。最后结果是[0 0 1 2 3 4 5 6 7 8]

对于只有 1 个单元格的滑动窗口来说似乎没问题，但它会让你大吃一惊：

>>> np.where((x[0,:]<5)&(x[0,:]>0),x[1,x[0,:]-2],x[1,:])
array([0, 9, 0, 1, 2, 5, 6, 7, 8, 9])

当您尝试通过两个单元格的窗口移动它时。

对于这个特定的问题，如果我们想把所有的东西都放在一行中，这样就可以了：

>>> for i in [1, 2, 3, 4, 5, 6]:
    print hstack((np.where(x[1,x[0,:]-i]<x[0, -i], x[1,x[0,:]-i], 0)[:5], x[0,5:]))

[0 0 1 2 3 5 6 7 8 9]
[0 0 0 1 2 5 6 7 8 9]
[0 0 0 0 1 5 6 7 8 9]
[0 0 0 0 0 5 6 7 8 9]
[0 0 0 0 0 5 6 7 8 9]
[0 0 0 0 0 5 6 7 8 9]

编辑： 现在我更好地理解了你原来的问题，基本上你想要一个二维数组并计算每个单元格周围的 N*N 单元格平均值。这很常见。首先，您可能希望将 N 限制为奇数，否则难以定义单元格周围的 2*2 平均值。假设我们想要 3*3 的平均值：

#In this example, the shape is (10,10)
>>> a1=\
array([[3, 7, 0, 9, 0, 8, 1, 4, 3, 3],
   [5, 6, 5, 2, 9, 2, 3, 5, 2, 9],
   [0, 9, 8, 5, 3, 1, 8, 1, 9, 4],
   [7, 4, 0, 0, 9, 3, 3, 3, 5, 4],
   [3, 1, 2, 4, 8, 8, 2, 1, 9, 6],
   [0, 0, 3, 9, 3, 0, 9, 1, 3, 3],
   [1, 2, 7, 4, 6, 6, 2, 6, 2, 1],
   [3, 9, 8, 5, 0, 3, 1, 4, 0, 5],
   [0, 3, 1, 4, 9, 9, 7, 5, 4, 5],
   [4, 3, 8, 7, 8, 6, 8, 1, 1, 8]])
#move your original array 'a1' around, use range(-2,2) for 5*5 average and so on
>>> movea1=[a1[np.clip(np.arange(10)+i, 0, 9)][:,np.clip(np.arange(10)+j, 0, 9)] for i, j in itertools.product(*[range(-1,2),]*2)]
#then just take the average
>>> averagea1=np.mean(np.array(movea1), axis=0)
#trim the result array, because the cells among the edges do not have 3*3 average
>>> averagea1[1:10-1, 1:10-1]
array([[ 4.77777778,  5.66666667,  4.55555556,  4.33333333,  3.88888889,
     3.66666667,  4.        ,  4.44444444],
   [ 4.88888889,  4.33333333,  4.55555556,  3.77777778,  4.55555556,
     3.22222222,  4.33333333,  4.66666667],
   [ 3.77777778,  3.66666667,  4.33333333,  4.55555556,  5.        ,
     3.33333333,  4.55555556,  4.66666667],
   [ 2.22222222,  2.55555556,  4.22222222,  4.88888889,  5.        ,
     3.33333333,  4.        ,  3.88888889],
   [ 2.11111111,  3.55555556,  5.11111111,  5.33333333,  4.88888889,
     3.88888889,  3.88888889,  3.55555556],
   [ 3.66666667,  5.22222222,  5.        ,  4.        ,  3.33333333,
     3.55555556,  3.11111111,  2.77777778],
   [ 3.77777778,  4.77777778,  4.88888889,  5.11111111,  4.77777778,
     4.77777778,  3.44444444,  3.55555556],
   [ 4.33333333,  5.33333333,  5.55555556,  5.66666667,  5.66666667,
     4.88888889,  3.44444444,  3.66666667]])

我认为您不需要将二维数组展平，这会导致混乱。此外，如果您想以不同方式处理边缘元素而不是仅仅将它们修剪掉，请考虑在“移动原始数组”步骤中使用 np.ma 制作掩码数组。