使用 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+ 加速！

利用numexpr 进行大数据计算以进一步提高性能。提升

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

所以，我们最终会得到这样的结果 -

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 了解如何控制多核功能。