python - 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,:])

索引 <5 的每个当前值的 n+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]]

有没有办法解决？我的方法完全错误吗？任何意见将不胜感激。

编辑 ：

这就是我想要实现的，我将一个矩阵展平为一个 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]的所有邻居的平均值。请注意，第一列和最后一列以及第一和最后一行受边界条件的约束，因此可能对您的应用程序无效（如果需要，您可以使用控制边界规则）。i,jarri,jmode=

因为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]])

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