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

Posted 2023-02-23

【中文标题】使用 strides 实现有效的移动平均滤波器

【英文标题】：Using strides for an efficient moving average filter

【发布时间】：2011-06-23 14:59:30

【问题描述】：

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

这是我目前所拥有的。它查看原始数组，然后将其滚动必要的数量，并对内核值求和以计算平均值。我知道边缘处理不正确，但我可以事后处理......有没有更好更快的方法？目标是过滤最大为 5000x5000 x 16 层的大型浮点数组，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)

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

当前代码：

使用 stride_tricks 生成一个类似 [[0,1,2],[1,2,3],[2,3,4]...] 的数组，它对应于过滤器内核的顶行。 沿纵轴滚动得到内核的中间行[[10,11,12],[11,12,13],[13,14,15]...]并将其添加到数组中我进去了 1) 重复获取内核的底部行 [[20,21,22],[21,22,23],[22,23,24]...]。在这一点上，我将每一行的总和除以过滤器中的元素数，得到每个像素的平均值，（移动 1 行和 1 列，边缘周围有一些奇怪的地方，但我可以稍后再处理）。

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

【问题讨论】：

我尝试运行您的代码，但出现内存损坏错误。我在 64 位 Ubuntu 10.10 上运行 Python 2.6.6 和 Numpy 1.3.0。错误看起来像*** glibc detected *** python: double free or corruption (!prev): 0x0000000002526d30 ***。

我能问一下为什么您使用浮点数（我假设是 64 位）来表示可以（可能）更有效地使用整数存储和计算的图像吗？

您的示例是 2D 数组，但您将数据描述为 3D。您是否对 16 层中的每一层都执行此操作？

@mtrw：我使用的是 Python 2.6.6。和 Windows XP SP3 上的 Numpy 1.4.1。我不知道那个错误是什么意思！

@Paul：数据是 3D（具有 16 个通道的图像），但可以作为单独的层过滤。我正在使用浮点数，因为这些值是雷达后向散射幅度值，并且不能选择截断或重新缩放。我最终将需要使用 Float32。

【参考方案1】：

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

@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/[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)

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

至少我的 0.02 美元...

【讨论】：

说得好：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....。当然，您对内存“爆炸”的解释很重要。您的回答（至少对我而言）的总结是：'不要太远钓鱼，quarenteed 渔获物已经在 scipy 中了'。谢谢

感谢乔，您的回答。在 rolling_window 中，if not hasattr(...): 是否应该返回 rolling_window_lastaxis(...) 而不是 rolling_window？

@unutbu - 完全正确！那是我的错别字...（我重命名了函数并忘记更改其中的那部分。）谢谢！

是否可以指定步长？

【参考方案2】：

我对 Python 不够熟悉，无法为此编写代码，但加速卷积的两种最佳方法是分离滤波器或使用傅里叶变换。

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

傅里叶变换：conv(A,B) 等价于ifft(fft(A)*fft(B))，即直接空间中的卷积变成傅里叶空间中的乘法，其中A 是您的图像，B 是您的过滤器.由于傅里叶变换的（逐元素）乘法要求 A 和 B 的大小相同，因此 B 是 size(A) 的数组，您的内核位于图像的最中心，其他位置为零。要将 3×3 内核放置在数组的中心，您可能必须将 A 填充为奇数大小。根据您对傅立叶变换的实现，这可能比卷积快很多（如果您多次应用相同的过滤器，您可以预先计算 fft(B)，从而节省 30% 的计算时间）。

【讨论】：

不管怎样，在python中，这些分别在scipy.ndimage.uniform_filter和scipy.signal.fftconvolve中实现。

@Jonas：酷！分离过滤器方法效果很好，正如您所说，随着内核大小的增加，它可以节省更多时间。对于 5000x5000 数组，在 11x11 内核大小下，使用 ndimage.convolve 的 2d 卷积得到 7.7 秒，使用 ndimage.convolve1d 得到两个 1d 卷积的 2.0 秒。对于您的第二个解决方案，B 是什么？

@Benjamin：我已经扩展了对第二种解决方案的解释

@Joe Kington：谢谢！如果我正确理解帮助，fftconvolve 不允许您预先计算 fft(B)，对吧？

@Benjamin - uniform_filter 已经做了重复的 1D 卷积，这是值得的。 @Jonas - 不，它没有......这只是一次说服的功能。

【参考方案3】：

让我们看看：

您的问题不是很清楚，但我现在假设您希望显着提高这种平均水平。

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)

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

更新： 首先，警告：当前状态的代码不能正确适应“内核”形状。然而，这不是我现在最关心的问题（无论如何，我已经有了如何正确适应的想法）。

我刚刚直观地选择了 4D A 的新形状，对我来说，考虑将 2D“内核”中心以原始 2D A 的每个网格位置为中心真的很有意义。

但 4D 塑形实际上可能不是“最好的”塑形。我认为这里真正的问题是求和的性能。应该能够找到（4D A 的）“最佳顺序”，以便充分利用您的机器缓存架构。但是，对于与您的计算机缓存“协作”的“小型”数组和不与计算机缓存“协作”的“小型”数组，该顺序可能不一样（至少不是那么简单的方式）。

更新 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)

【讨论】：

@eat：这很有趣。我看到您已经解决了我的边缘问题，尽管您的过滤器大小是硬编码的；）。在您的 as_strided 行中，应该是 n, m 而不是 n, n？

@eat：我稍微修改了你的代码，它运行良好。不过，我无法理解正在发生的事情。您能否描述一下 as_strided 线在做什么以及为什么要选择这些形状和步幅值？

@eat: 这些 n 中的一个不应该是 m 吗？

@Bejamin，@Paul：是的，有些错别字，只是在我的答案中更改了这些。无论如何，我认为这里真正的烫手山芋是：我们可以期待这个实现有显着的改进吗？稍后我将尝试进一步澄清我的答案。谢谢

@eat：即使在 3x3 内核大小的情况下，似乎也存在较大数组（如 5000x5000）的问题...

【参考方案4】：

我确信需要修复的一件事是您的视图数组b。

它有一些来自未分配内存的项目，所以你会遇到崩溃。

鉴于您对算法的新描述，需要修复的第一件事是您超出了a 的分配范围：

bshape = (a.size-filtsize+1, filtsize)
bstrides = (a.itemsize, a.itemsize)
b = numpy.lib.stride_tricks.as_strided(a, shape=bshape, strides=bstrides)

---

更新

因为我还没有完全掌握方法，而且似乎有更简单的方法可以解决问题，所以我只是把它放在这里：

A = numpy.arange(100).reshape((10,10))

shifts = [(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]
B = A[1:-1, 1:-1].copy()
for dx,dy in shifts:
    xstop = -1+dx or None
    ystop = -1+dy or None
    B += A[1+dx:xstop, 1+dy:ystop]
B /= 9

...这似乎是一种直截了当的方法。唯一无关的操作是它只分配和填充B 一次。无论如何，所有的加法、除法和索引都必须完成。如果你正在做 16 个波段，如果你的意图是保存图像，你仍然只需要分配一次B。即使这没有帮助，它也可以澄清为什么我不理解这个问题，或者至少可以作为其他方法加速时间的基准。这在我的笔记本电脑上的 5k x 5k 数组 float64 上运行 2.6 秒，其中 0.5 是 B 的创建

【讨论】：

我刚刚计时了，你的方法比我的快 10 倍。该比率似乎相当恒定（即它不取决于输入大小）。一个仓促的结论是 stride_tricks 有时是有用的技巧，但它们不需要任何性能提升？ Alltough 可能存在一些其他技巧，而不是我的表现更好。谢谢

@eat：对于 5000x5000 数组和 3x3 过滤器，我将 @eat 的结果计时 3.9 秒，@Paul 的结果计时 1.9 秒，scipy.ndimage.filters.convolve 计时 1.4 秒。在该数组大小下，strides 解决方案不适用于较大的内核大小。我将升级@Paul 的解决方案以接受可变内核大小并进行比较。但似乎 scipy.ndimage.filters.convolve 仍然是最快的解决方案......

如何处理二维数组中的 NaN 值？