使用步幅获得有效的移动平均滤波器
image-processing
numpy
python
5
0

我最近在这篇帖子答案中了解了步幅 ,并且想知道如何使用它们来计算比我在本文中建议 (使用卷积滤波器)更有效的移动平均滤波器。

到目前为止,这就是我所拥有的。它查看原始数组,然后将其滚动所需的量,并对内核值求和以计算平均值。我知道边缘处理不正确,但是以后我可以解决...是否有更好,更快的方法?目的是过滤大小最大为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个值或内核元素的总和,或者有人可以说服我使用另一种更有效的方法...

参考资料:
Stack Overflow
收藏
评论
共 4 个回答
高赞 时间 活跃

对于它的价值,这是您将使用“花式”大踏步花样来做的事情。我本打算昨天发布这个,但由于实际工作而分心! :)

@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维数组的选定轴进行操作(即,在4维数组上的filtsize = (0,3,0,3)可以使我们得到6维视图)。

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

但是,由于在mean (或std等)的每一步上存储的临时数组都比原始数组大得多,因此这根本没有提高内存效率!同样,它也不会很快。

ndimage的等效ndimage只是:

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

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

无论如何,只有我的$ 0.02 ...

收藏
评论

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

分离的滤镜 :卷积为O(M * N),其中M和N分别是图像和滤镜中的像素数。由于使用3-by-3内核进行平均过滤等效于首先使用3-by-1内核进行过滤,然后再使用1-by-3内核进行过滤,因此可以得到(3+3)/(3*3) =〜通过使用两个1-d内核进行连续卷积,速度提高了30%(随着内核变大,这显然会变得更好)。当然,您可能仍然可以在这里使用大步技巧。

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

收藏
评论

让我们来看看:

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

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)
收藏
评论

我确信需要解决的一件事是您的视图数组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即使这无济于事,也可以弄清楚为什么我不理解这个问题,或者至少可以作为衡量其他方法加速的基准。这在我的笔记本电脑上以2.6秒的速度在5k x 5k的float64数组上运行,其中0.5是B的创建

收藏
评论
新手导航
  • 社区规范
  • 提出问题
  • 进行投票
  • 个人资料
  • 优化问题
  • 回答问题

关于我们

常见问题

内容许可

联系我们

@2020 AskGo
京ICP备20001863号