提高阵列范围过滤器中的内存使用率,以避免块处理

Question

我正在实现一些卫星图像滤波器,从一个称为增强型Lee滤波器开始.图像很容易达到5000x5000像素甚至更多.我当前的实现是内存耗尽,试图计算那些大型数组上的过滤器(请注意,移动平均值和移动stddev过滤器可以一次运行).主要困难是为了返回最终过滤的数组,必须保留在内存中的数组的数量.在这个问题中,我请求帮助改进块处理函数,但我的问题是:有没有办法改进这段代码,这样我就不需要使用块处理了？

    def moving_average(Ic, filtsize):
        Im = numpy.empty(Ic.shape, dtype='Float32')
        scipy.ndimage.filters.uniform_filter(Ic, filtsize, output=Im)
        return Im

    def moving_stddev(Ic, filtsize):
        Im = numpy.empty(Ic.shape, dtype='Float32')
        scipy.ndimage.filters.uniform_filter(Ic, filtsize, output=Im)
        S = numpy.empty(Ic.shape, dtype='Float32')
        scipy.ndimage.filters.uniform_filter(((Ic-Im) ** 2), filtsize, output=S)
        return numpy.sqrt(S)

    def enh_lee(Ic, filtsize, nlooks, dfactor):
        # Implementation based on PCI Geomatica's FELEE function documentation
        Ci = moving_stddev(Ic, filtsize) / moving_average(Ic, filtsize) #1st array in memory
        Cu = numpy.sqrt(1 / nlooks) #scalar
        Cmax = numpy.sqrt(1 + (2 * nlooks)) #scalar
        W = numpy.exp(-dfactor * (Ci - Cu) / (Cmax - Ci)) #2nd array in memory
        Im = moving_average(Ic, filtsize) #3rd array in memory
        If = Im * W + Ic * (1 - W) #4th array in memory
        W = None # Back to 3 arrays in memory
        return numpy.select([Ci <= Cu, (Cu < Ci) * (Ci < Cmax), Ci >= Cmax], [Im, If, Ic])

where nlooks和dfactorscalars Ic是未过滤的数组.

根据你的建议进行编辑(我也在查看numexpr),我的改进代码enh_lee如下,但仍然不足以超越最后一步而不会耗尽内存:

def enh_lee(Ic, filtsize, nlooks, dfactor):
    Im = moving_average(Ic, filtsize)
    Ci = moving_stddev(Ic, filtsize)
    Ci /= Im
    Cu = numpy.sqrt(1 / nlooks)
    Cmax = numpy.sqrt(1 + (2 * nlooks))

    W = Ci
    W -= Cu
    W /= Cmax - Ci
    W *= -dfactor
    numpy.exp(W, W)

    If = 1
    If -= W
    If *= Ic
    If += Im * W
    W = None
    return numpy.select([Ci <= Cu, (Cu < Ci) & (Ci < Cmax), Ci >= Cmax], [Im, If, Ic])

Answer 1

你可以在这里做几个(内存使用)优化...要记住的一些技巧是:

1. 大多数numpy函数都使用一个out参数,而不是用于指定输出数组而不是返回副本.例如,np.sqrt(x, x)将采用数组的平方根.
2. x += 1使用一半的内存x = x + 1,因为后者制作临时副本.如果有可能,尽量向上拆分成计算*=,+=,/=,等.当它不是,使用numexpr,作为@eumiro建议.(或者只是使用了numexpr ...在许多情况下它非常方便.)

所以,首先,这是你的原始函数的性能,10000x10000随机数据阵列和filtsize3:

原始功能的内存使用情况

有趣的是最后的大峰值.这些都发生在你的numpy.select(...)位.有很多地方你无意中创建了额外的临时阵列,但它们大多无关紧要,因为它们在通话过程中不堪重负select.

无论如何,如果我们用这个相当冗长的版本替换你的原始(干净和紧凑)代码,你可以显着优化你的内存使用:

import numpy
import scipy.ndimage

def main(x=None):
    if x is None:
        ni, nj = 10000, 10000
        x = numpy.arange(ni*nj, dtype=numpy.float32).reshape(ni,nj)
    filtsize = 3
    nlooks = 10.0
    dfactor = 10.0
    x = enh_lee(x, filtsize, nlooks, dfactor)
    return x

def moving_average(Ic, filtsize):
    Im = numpy.empty(Ic.shape, dtype='Float32')
    scipy.ndimage.filters.uniform_filter(Ic, filtsize, output=Im)
    return Im

def moving_stddev(Ic, filtsize):
    Im = numpy.empty(Ic.shape, dtype='Float32')
    scipy.ndimage.filters.uniform_filter(Ic, filtsize, output=Im)
    Im *= -1
    Im += Ic
    Im **= 2
    scipy.ndimage.filters.uniform_filter(Im, filtsize, output=Im)
    return numpy.sqrt(Im, Im)

def enh_lee(Ic, filtsize, nlooks, dfactor):
    # Implementation based on PCI Geomatica's FELEE function documentation
    Ci = moving_stddev(Ic, filtsize) 
    Im = moving_average(Ic, filtsize) 
    Ci /= Im

    Cu = numpy.sqrt(1 / nlooks).astype(numpy.float32) 
    Cmax = numpy.sqrt(1 + (2 * nlooks)).astype(numpy.float32) 

    W = Ci.copy()
    W -= Cu
    W *= -dfactor
    W /= Cmax - Ci
    W = numpy.exp(W, W)

    If = Im * W
    W *= -1
    W += 1
    W *= Ic
    If += W
    del W

    # Replace the call to numpy.select
    out = If
    filter = Ci <= Cu
    numpy.putmask(out, filter, Im)
    del Im

    filter = Ci >= Cmax
    numpy.putmask(out, filter, Ic)
    return out

if __name__ == '__main__':
    main()

这是此代码的结果内存配置文件:

基于Numpy的优化版本的内存使用情况

因此,我们大大减少了内存使用量,但代码的可读性稍差(imo).

然而,最后三个峰值是两个numpy.where电话......

如果numpy.where采用out参数,我们可以进一步将峰值内存使用量减少大约300Mb左右.不幸的是,它没有,我不知道更有效的内存方式...

我们可以使用numpy.putmask替换呼叫numpy.select并就地进行操作(感谢@eumiro 在一个完全不同的问题中提到这一点.)

如果我们使用numexpr优化事物,我们会得到相当清晰的代码(与上面的纯粹numpy版本相比,而不是原始代码).在这个版本中你可能会减少一点内存使用量...我不是非常熟悉numexpr,除了使用它几次之外.

import numpy
import scipy.ndimage
import numexpr as ne

def main(x=None):
    if x is None:
        ni, nj = 10000, 10000
        x = numpy.arange(ni*nj, dtype=numpy.float32).reshape(ni,nj)
    filtsize = 3
    nlooks = 10.0
    dfactor = 10.0
    x = enh_lee(x, filtsize, nlooks, dfactor)
    return x

def moving_average(Ic, filtsize):
    Im = numpy.empty(Ic.shape, dtype='Float32')
    scipy.ndimage.filters.uniform_filter(Ic, filtsize, output=Im)
    return Im

def moving_stddev(Ic, filtsize):
    Im = numpy.empty(Ic.shape, dtype='Float32')
    scipy.ndimage.filters.uniform_filter(Ic, filtsize, output=Im)
    Im = ne.evaluate('((Ic-Im) ** 2)')
    scipy.ndimage.filters.uniform_filter(Im, filtsize, output=Im)
    return ne.evaluate('sqrt(Im)')

def enh_lee(Ic, filtsize, nlooks, dfactor):
    # Implementation based on PCI Geomatica's FELEE function documentation
    Ci = moving_stddev(Ic, filtsize) 
    Im = moving_average(Ic, filtsize) 
    Ci /= Im

    Cu = numpy.sqrt(1 / nlooks).astype(numpy.float32) 
    Cmax = numpy.sqrt(1 + (2 * nlooks)).astype(numpy.float32) 

    W = ne.evaluate('exp(-dfactor * (Ci - Cu) / (Cmax - Ci))') 
    If = ne.evaluate('Im * W + Ic * (1 - W)') 
    del W

    out = ne.evaluate('where(Ci <= Cu, Im, If)') 
    del Im
    del If

    out = ne.evaluate('where(Ci >= Cmax, Ic, out)')
    return out

if __name__ == '__main__':
    main()

这里是numexpr版本的内存使用情况:(请注意,与原始版本相比,执行时间减少了一半!)

基于Numexpr的优化版本的内存使用情况*

在调用where(替换调用select)期间仍然使用最大的内存.但是,峰值内存使用量已大幅减少.进一步减少它的最简单方法是找到一些方法select在其中一个阵列上就地操作.使用cython执行此操作相当容易(嵌套循环在纯python中会相当慢,并且numpy中的任何类型的布尔索引都会创建一个额外的副本).通过简单地按照你一直在做的方式对输入数组进行分块,你可能会更好,不过......

与旁注一样,更新版本产生的输出与原始代码相同.原始的基于numpy的代码中有一个拼写错误...