利用步幅实现有效的移动平均fi

Question

我最近了解了answer to this post中的strides，并想知道如何使用它们比我在this post（使用卷积滤波器）中建议的更有效地计算移动平均滤波器。

这就是我目前所拥有的。它获取原始数组的视图，然后将其按所需数量滚动，并对内核值求和以计算平均值。我知道边缘处理不正确，但我可以在以后处理。。。有更好更快的方法吗？其目标是过滤大小高达5000x5000 x 16层的大型浮点数组，这项任务的速度相当慢。

请注意，我正在寻找8邻连接性，即3x3滤波器取9个像素的平均值（焦点像素周围8个像素），并将该值分配给新图像中的像素。

import numpy, scipy

filtsize = 3
a = numpy.arange(100).reshape((10,10))
b = numpy.lib.stride_tricks.as_strided(a, shape=(a.size,filtsize), strides=(a.itemsize, a.itemsize))
for i in range(0, filtsize-1):
    if i > 0:
        b += numpy.roll(b, -(pow(filtsize,2)+1)*i, 0)
filtered = (numpy.sum(b, 1) / pow(filtsize,2)).reshape((a.shape[0],a.shape[1]))
scipy.misc.imsave("average.jpg", filtered)

编辑关于我如何看待此工作的说明：

当前代码：

1. 使用stride_技巧生成一个数组，如[[0,1,2]，[1,2,3]，[2,3,4]…]，它对应于过滤器内核的顶行。
2. 沿垂直轴滚动以获得内核的中间行[[10,11,12]，[11,12,13]，[13,14,15]…]并将其添加到我在1中获得的数组中）
3. 重复以获取内核的最下面一行[[20,21,22]，[21,22,23]，[22,23,24]…]。在这一点上，我取每一行的和除以过滤器中的元素数，得到每个像素的平均值（移动1行和1列，边缘有一些奇怪的地方，但我可以稍后处理）。

我所希望的是更好地使用stride_技巧，直接获取整个数组的9个值或内核元素的总和，或者有人能说服我使用另一种更有效的方法。。。

Answer 1

让我们看看：

你的问题还不太清楚，但我现在假设你想大大提高这种平均数。

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

def mf(A, k_shape= (3, 3)):
    m= A.shape[0]- 2
    n= A.shape[1]- 2
    strides= A.strides+ A.strides
    new_shape= (m, n, k_shape[0], k_shape[1])
    A= st.as_strided(A, shape= new_shape, strides= strides)
    return np.sum(np.sum(A, -1), -1)/ np.prod(k_shape)

if __name__ == '__main__':
    A= np.arange(100).reshape((10, 10))
    print mf(A)

现在，您实际期望的性能改进是什么？

更新：
首先，有一个警告：处于当前状态的代码不能正确地适应“内核”形状。不过，这并不是我现在最关心的问题（不管怎样，我的想法是已经准备好如何适当地适应了）。

我刚刚直观地选择了一个4da的新形状，对我来说，考虑一个2D‘核’中心以原始2da的每个网格位置为中心是非常有意义的

但4D造型可能并不是最好的。我认为这里真正的问题是求和的表现。我们应该能够找到“最佳顺序”（4da），以便充分利用您的机器缓存体系结构。但是，对于“小”数组（这种数组与机器缓存“协同工作”）和那些“大”数组（至少不是那么直接的方式）来说，顺序可能不一样。

更新2:
这是mf的一个稍微修改的版本。很明显，最好先重塑为一个3D数组，然后不求和，只做点积（这样做的好处是，内核可以是任意的）。但是它仍然比Pauls更新的功能慢3倍（在我的机器上）。

def mf(A):
    k_shape= (3, 3)
    k= np.prod(k_shape)
    m= A.shape[0]- 2
    n= A.shape[1]- 2
    strides= A.strides* 2
    new_shape= (m, n)+ k_shape
    A= st.as_strided(A, shape= new_shape, strides= strides)
    w= np.ones(k)/ k
    return np.dot(A.reshape((m, n, -1)), w)