利用步幅实现有效的移动平均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)

Answer 2

我对Python还不太熟悉，因此无法编写代码，但加速卷积的两个最佳方法是分离滤波器或使用Fourier变换。

分离滤波器：卷积是O（M*N），其中M和N分别是图像和滤波器中的像素数。由于使用3×3内核的平均过滤相当于先使用3×1内核过滤，然后使用1×3内核过滤，因此通过使用两个1-d内核的连续卷积可以获得(3+3)/(3*3)=~30%的速度提高（这显然随着内核变大而变得更好）。当然，在这里你仍然可以使用跨步技巧。

Fourier变换等价于ifft(fft(A)*fft(B))，即直接空间中的卷积变成Fourier空间中的乘法，其中A是图像，B是滤波器。由于Fourier变换的（按元素）乘法要求A和B的大小相同，所以B是一个size(A)数组，内核位于图像的中心，其他地方都是零。要将3×3内核放在数组的中心，可能需要将A填充到奇数大小。根据Fourier变换的实现，这可能比卷积快得多（如果多次应用相同的滤波器，则可以预先计算fft(B)，再节省30%的计算时间）。

Answer 3

值得一提的是，以下是你如何使用“花式”跨步技巧。我昨天本来要发这个的，但是被实际工作搞得心烦意乱！：）

@Paul&；@eat都有很好的实现，使用各种其他方法来实现。为了继续前面问题的内容，我想我应该发布N维等价物。

但是，对于>；1D数组，您将无法显著优于scipy.ndimage函数。（scipy.ndimage.uniform_filter应该优于scipy.ndimage.convolve）

此外，如果您试图获得多维移动窗口，则在无意中复制数组时，可能会导致内存使用量爆炸。虽然最初的“滚动”数组只是原始数组内存中的一个视图，但是复制该数组的任何中间步骤都将生成一个比原始数组大个数量级的副本（例如，假设您使用的是100x100原始数组。。。它的视图（对于（3,3）大小的过滤器）将是98x98x3x3，但使用与原始内存相同的内存。但是，任何副本都将使用完全98x98x3x3数组所需的内存量！！）

基本上，当你想在ndarray的单轴上对移动窗口操作进行矢量化时，使用疯狂的跨步技巧是很好的。它使得计算移动标准差等事情变得非常容易，而且开销很小。当你想沿着多个轴开始这样做的时候，这是可能的，但是你通常最好使用更专门的函数。（如scipy.ndimage等）

无论如何，以下是你的做法：

import numpy as np

def rolling_window_lastaxis(a, window):
    """Directly taken from Erik Rigtorp's post to numpy-discussion.
    <http://www.mail-archive.com/numpy-discussion@scipy.org/msg29450.html>"""
    if window < 1:
       raise ValueError, "`window` must be at least 1."
    if window > a.shape[-1]:
       raise ValueError, "`window` is too long."
    shape = a.shape[:-1] + (a.shape[-1] - window + 1, window)
    strides = a.strides + (a.strides[-1],)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

def rolling_window(a, window):
    if not hasattr(window, '__iter__'):
        return rolling_window_lastaxis(a, window)
    for i, win in enumerate(window):
        if win > 1:
            a = a.swapaxes(i, -1)
            a = rolling_window_lastaxis(a, win)
            a = a.swapaxes(-2, i)
    return a

filtsize = (3, 3)
a = np.zeros((10,10), dtype=np.float)
a[5:7,5] = 1

b = rolling_window(a, filtsize)
blurred = b.mean(axis=-1).mean(axis=-1)

所以当我们做b = rolling_window(a, filtsize)的时候，我们得到的是一个8x8x3x3数组，它实际上是一个进入原始10x10数组相同内存的视图。我们可以很容易地沿着不同的轴使用不同的滤波器大小，或者只沿着N维数组的选定轴操作（即，在四维数组上的filtsize = (0,3,0,3)将为我们提供一个6维视图）。

然后，我们可以对最后一个轴重复应用任意函数，以有效地计算移动窗口中的内容。

但是，因为我们在mean（或std或其他任何步骤）的每一步都存储比原始数组大得多的临时数组，所以这根本不节省内存！也不会太快。

与ndimage等价的是：

blurred = scipy.ndimage.uniform_filter(a, filtsize, output=a)

这将处理各种边界条件，在不需要阵列临时副本的情况下进行“模糊”处理，并且可以非常快地veryfast。跨步技巧是将函数应用于沿一个轴移动的窗口的好方法，但它们不是沿多个轴移动的好方法，通常。。。。

不管怎样，只要我的0.02美元。。。