python - 使用 strides 实现有效的移动平均滤波器

Question

我最近了解到 strides在 answer to this post ，并且想知道如何使用它们比我在 this post 中提出的更有效地计算移动平均滤波器(使用卷积过滤器)。

这是我目前所拥有的。它查看原始数组，然后将其滚动必要的数量，并对内核值求和以计算平均值。我知道边缘处理不正确，但我可以事后处理......有没有更好更快的方法？目标是过滤大小高达 5000x5000 x 16 层的大型 float 组，scipy.ndimage.filters.convolve 的任务相当慢。

请注意，我正在寻找 8-neighbour 连接，即 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 列，边缘周围有一些奇怪的地方，但我可以稍后再处理)。

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

Answer 1

对于它的值(value)，这里是你如何使用“花哨的”跨步技巧来做到这一点。昨天本来打算发这个的，但是被实际工作分心了! :)

@Paul 和 @eat 都有很好的实现，使用各种其他方式来实现。只是为了继续前面的问题，我想我会发布 N 维等价物。

但是，对于 >1D 数组，您将无法显着击败 scipy.ndimage 函数。 (虽然 scipy.ndimage.uniform_filter 应该击败 scipy.ndimage.convolve)

此外，如果您试图获得一个多维移动窗口，那么每当您无意中复制了数组时，就有可能导致内存使用量激增。虽然最初的“滚动”数组只是原始数组内存的一个 View ，但复制数组的任何中间步骤都会生成一个比原始数组大数量级的副本(即让我们假设您正在使用一个 100x100 的原始数组...对其的 View (对于 (3,3) 的过滤器大小)将是 98x98x3x3，但使用与原始数组相同的内存。但是，任何副本都将使用该数量full 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 数组相同内存的 View 。我们可以很容易地在不同的轴上使用不同的过滤器大小，或者只在 N 维数组的选定轴上操作(即 filtsize = (0,3,0,3) 在 4 维数组会给我们一个 6 维 View )。

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

但是，因为我们在 mean(或 std 或其他任何步骤)的每个步骤中存储的临时数组都比原始数组大得多，所以这不是所有内存高效!它也不会很快。

ndimage 的等价物只是:

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

这将处理各种边界条件，就地“模糊”而不需要数组的临时副本，并且非常快。跨步技巧是沿 one 轴将函数应用于移动窗口的好方法，但它们不是沿多个轴执行此操作的好方法，通常......

至少我的 0.02 美元...