大步迈向高效的移动平均滤波器

Question

我最近在这篇文章的答案中了解了步幅，并且想知道如何使用它们比我在这篇文章中提出的)（使用卷积滤波器）更有效地计算移动平均滤波器。

到目前为止，这就是我所拥有的。它查看原始数组，然后将其滚动所需的量，然后对内核值求和以计算平均值。我知道边缘处理不正确，但是以后我可以解决…是否有更好，更快的方法？目标是过滤大小最大为5000x5000
x 16层的大型浮点数组，这scipy.ndimage.filters.convolve是一个相当慢的任务。

请注意，我正在寻找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_tricks生成类似于[[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 col，边缘周围有些奇怪，但我可以稍后再处理）。

我希望可以更好地利用stride_tricks来直接获取整个数组的9个值或内核元素的总和，或者有人可以说服我采用另一种更有效的方法…

Answer 1

对于它的价值，这是您将使用“花哨”的大步小技巧来做的事情。我本打算昨天发布这个，但是由于实际工作而分心！:)

@Paul和@eat都使用各种其他方法来实现良好的实现。为了继续前面的问题，我以为我会发布N维等效项。

但是，您将无法显着击败scipy.ndimage1D阵列的功能。（虽然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/[e-mail]/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)，在4维数组上将为我们提供6维视图）。

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

但是，由于我们在mean（std或其他任何步骤）的每个步骤中都存储了比原始数组大得多的临时数组，因此这根本没有提高内存效率！同样，它也不会很快。

等效于ndimage：

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

这将处理各种边界条件，在不需要阵列的临时副本的情况下就地进行“模糊处理”，而且速度 非常 快。跨步技巧是将功能沿 一个
轴应用于移动窗口的一种好方法，但通常不是沿多个轴进行操作的一种好方法。

无论如何，只要$ 0.02美元…