使用Numpy可以更有效地生成过滤器

Question

我正在生成一些过滤器，用于实现快速傅里叶变换模糊和锐化图像操作。正确生成过滤器，但计算持续时间很长。

我当前生成过滤器的方式是逐项迭代所需过滤器的尺寸。我明白我需要使用Numpy来解决这个问题，但我无法弄清楚究竟是怎么回事。这是我生成高斯滤波器的代码：

def gaussian_filter(mode, size, cutoff):
    filterImage = np.zeros(size, np.float64)    
    cutoffTerm = 2 * (cutoff ** 2)
    v = np.asarray([size[0] // 2, size[1] // 2])

    for px in range(0, size[0]):
        for py in range(0, size[1]):
            u = np.asarray([px, py])
            Duv = np.linalg.norm(u - v)
            distance = -1 * (Duv ** 2)
            result = pow(np.e, distance / cutoffTerm)
            if mode == 'low':
                filterImage.itemset((px, py), result)
            elif mode == 'high':
                filterImage.itemset((px, py), 1 - result)

    return filterImage

生成大小为1920 x 1080的过滤器需要70.36秒，这是完全不可接受的。任何想法将不胜感激。

Answer 1

这是一个利用broadcasting的矢量化的 -

def gaussian_filter_vectorized(mode, size, cutoff):
    cutoffTerm = 2 * (cutoff ** 2)
    v = np.asarray([size[0] // 2, size[1] // 2])

    I,J = np.ogrid[:size[0],:size[1]]
    p,q = I-v[0],J-v[1]
    Dsq = p**2 + q**2
    d = -1 * Dsq
    R = np.power(np.e,d/cutoffTerm)
    if mode == 'low':
        return R
    elif mode == 'high':
        return 1-R

关于大size的定时 -

In [80]: N = 100
    ...: %timeit gaussian_filter(mode='low', size=(N,N), cutoff=N)
    ...: %timeit gaussian_filter_vectorized(mode='low', size=(N,N), cutoff=N)
10 loops, best of 3: 65.2 ms per loop
1000 loops, best of 3: 225 µs per loop

In [81]: N = 1000
    ...: %timeit gaussian_filter(mode='low', size=(N,N), cutoff=N)
    ...: %timeit gaussian_filter_vectorized(mode='low', size=(N,N), cutoff=N)
1 loop, best of 3: 6.5 s per loop
10 loops, best of 3: 29.8 ms per loop

200x+加速！

Leverage numexpr on large data computations for further perf. boost

在处理大数据时，如果预期的操作可以表示为算术运算，我们也可以使用支持多核处理的numexpr module。为了解决我们的情况，我们可以使用Dsq = p**2 + q**2函数替换R = np.power(np.e,d/cutoffTerm)和numexpr中的步骤与numexpr.evaluate等效的步骤。

所以，我们最终会得到这样的东西 -

import numexpr as ne

def gaussian_filter_vectorized_numexpr(mode, size, cutoff):
    cutoffTerm = 2 * (cutoff ** 2)

    I,J = np.ogrid[:size[0],:size[1]]
    v0,v1 = size[0] // 2, size[1] // 2
    p,q = I-v0,J-v1    
    E = np.e
    if mode == 'low':
        return ne.evaluate('E**(-1*(p**2+q**2)/cutoffTerm)')
    elif mode == 'high':
        return ne.evaluate('1-E**(-1*(p**2+q**2)/cutoffTerm)')

关于1920x1080尺寸图像的计时 -

In [2]: M,N=1920,1080
   ...: %timeit gaussian_filter(mode='low', size=(M,N), cutoff=N)
   ...: %timeit gaussian_filter_vectorized(mode='low', size=(M,N), cutoff=N)
   ...: %timeit gaussian_filter_vectorized_numexpr(mode='low',size=(M,N),cutoff=N)
1 loop, best of 3: 13.9 s per loop
10 loops, best of 3: 63.3 ms per loop
100 loops, best of 3: 9.48 ms per loop

接近1500x加速在这里！

这是与8线程。因此，随着可用于计算的更多线程，它应该进一步改进。 Related post关于如何控制多核功能。