使用 strides 实现高效的移动平均滤波器 [英] Using strides for an efficient moving average filter

问题描述

我最近了解了 strides这篇文章的答案，我想知道我怎么能使用它们比我在 这篇文章(使用卷积滤波器).

I recently learned about strides in the answer to this post, and was wondering how I could use them to compute a moving average filter more efficiently than what I proposed in this post (using convolution filters).

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

This is what I have so far. It takes a view of the original array then rolls it by the necessary amount and sums the kernel values to compute the average. I am aware that the edges are not handled correctly, but I can take care of that afterward... Is there a better and faster way? The objective is to filter large floating point arrays up to 5000x5000 x 16 layers in size, a task that scipy.ndimage.filters.convolve is fairly slow at.

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

Note that I am looking for 8-neighbour connectivity, that is a 3x3 filter takes the average of 9 pixels (8 around the focal pixel) and assigns that value to the pixel in the new image.

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)

编辑澄清我如何看待这个工作:

EDIT Clarification on how I see this working:

当前代码:

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 个值或内核元素的总和，或者有人可以说服我采用另一种更有效的方法...

What I was hoping for is a better use of stride_tricks to get the 9 values or the sum of the kernel elements directly, for the entire array, or that someone can convince me of another more efficient method...

推荐答案

就其价值而言，以下是您如何使用花哨的"跨步技巧来做到这一点.我昨天打算发这个，但被实际工作分散了注意力！:)

For what it's worth, here's how you'd do it using "fancy" striding tricks. I was going to post this yesterday, but got distracted by actual work! :)

@保罗&@eat 两者都使用其他各种方式实现了很好的实现.只是为了继续之前的问题，我想我会发布 N 维等价物.

@Paul & @eat both have nice implementations using various other ways of doing this. Just to continue things from the earlier question, I figured I'd post the N-dimensional equivalent.

但是，对于 >1D 数组，您将无法明显击败 scipy.ndimage 函数.(不过，scipy.ndimage.uniform_filter 应该击败 scipy.ndimage.convolve)

You're not going to be able to significantly beat scipy.ndimage functions for >1D arrays, however. (scipy.ndimage.uniform_filter should beat scipy.ndimage.convolve, though)

此外，如果您想获得一个多维移动窗口，那么每当您不经意地复制数组时，就有可能导致内存使用量激增.虽然最初的滚动"数组只是原始数组内存中的一个视图，但复制数组的任何中间步骤都会制作一个比原始数组大数量级的副本(即让我们假设你正在使用一个 100x100 的原始数组......进入它的视图(对于 (3,3) 的过滤器大小)将是 98x98x3x3 但使用与原始数组相同的内存.但是，任何副本都将使用该数量一个完整 98x98x3x3 数组所需的内存！！)

Moreover, if you're trying to get a multidimensional moving window, you risk having memory usage blow up whenever you inadvertently make a copy of your array. While the initial "rolling" array is just a view into the memory of your original array, any intermediate steps that copy the array will make a copy that is orders of magnitude larger than your original array (i.e. Let's say that you're working with a 100x100 original array... The view into it (for a filter size of (3,3)) will be 98x98x3x3 but use the same memory as the original. However, any copies will use the amount of memory that a full 98x98x3x3 array would!!)

基本上，当您想在 ndarray 的单轴上矢量化移动窗口操作时，使用疯狂的跨步技巧非常有用.它使得计算诸如移动标准偏差之类的东西变得非常容易，而且开销很小.当您想开始沿多个轴执行此操作时，这是可能的，但通常最好使用更专业的函数.(如scipy.ndimage等)

Basically, using crazy striding tricks is great for when you want to vectorize moving window operations on a single axis of an ndarray. It makes it really easy to calculate things like a moving standard deviation, etc with very little overhead. When you want to start doing this along multiple axes, it's possible, but you're usually better off with more specialized functions. (Such as scipy.ndimage, etc)

无论如何，这就是你的做法:

At any rate, here's how you do it:

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 维的视图).

So what we get when we do b = rolling_window(a, filtsize) is an 8x8x3x3 array, that's actually a view into the same memory as the original 10x10 array. We could have just as easily used different filter size along different axes or operated only along selected axes of an N-dimensional array (i.e. filtsize = (0,3,0,3) on a 4-dimensional array would give us a 6 dimensional view).

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

We can then apply an arbitrary function to the last axis repeatedly to effectively calculate things in a moving window.

然而，因为我们在 mean(或 std 或其他)的每个步骤中存储的临时数组比原始数组大得多，所以这不是所有的内存效率！它也不会非常快.

However, because we're storing temporary arrays that are much bigger than our original array on each step of mean (or std or whatever), this is not at all memory efficient! It's also not going to be terribly fast, either.

ndimage 的等价物只是:

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

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

This will handle a variety of boundary conditions, do the "blurring" in-place without requiring a temporary copy of the array, and be very fast. Striding tricks are a good way to apply a function to a moving window along one axis, but they're not a good way to do it along multiple axes, usually....

无论如何，只要我的 0.02 美元......

Just my $0.02, at any rate...

这篇关于使用 strides 实现高效的移动平均滤波器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！