python - 改进阵列范围过滤器中的内存使用以避免 block 处理

Question

我正在实现一些卫星图像过滤器，首先是一种称为 Enhanced Lee 过滤器的过滤器。图像很容易达到 5000x5000 像素甚至更多。我当前的实现在尝试计算那些大型数组上的过滤器时内存不足(请注意，移动平均和移动标准开发过滤器可以一次性运行)。主要困难是为了返回最终过滤后的数组必须保存在内存中的数组数量。在 this question ，我寻求有关改进 block 处理功能的帮助，但我的问题是:是否有改进此代码以便我不需要使用 block 处理的方法？

    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])

其中 nlooks 和 dfactor 是标量，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 随机数据数组和 3 的 filtsize 时原始函数的性能:

原始函数的内存使用情况

有趣的是末端的大尖峰。这些发生在您的 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的优化版内存使用情况

因此，我们大大减少了内存使用量，但代码的可读性稍差(i.m.o.)。

但是，最后三个峰值是两个 numpy.where 调用...

如果 numpy.where 有一个 out 参数，我们可以进一步将峰值内存使用量再减少 ~300Mb 左右。不幸的是，它没有，而且我不知道有更节省内存的方法......

我们可以使用 numpy.putmask 来替换对 numpy.select 的调用，并就地进行操作(感谢@eumiro mentioning this in an entirely different question。)

如果我们使用 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 中的任何类型的 bool 索引都会创建一个附加副本)。不过，像您一直在做的那样简单地对输入数组进行分块可能会更好……

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