Python中最快的2D卷积或图像滤波器

Question

有很多用户在问关于在numpy或scipy中进行图像卷积时的速度或内存消耗的问题 [1, 2, 3, 4]。根据回复和我使用Numpy的经验，我觉得这可能是numpy相比于Matlab或IDL的一个主要短板。

到目前为止，没有任何答案真正解决了整体问题，所以我想问：“在Python中计算2D卷积的最快方法是什么？”常见的Python模块都可以使用：numpy、scipy和PIL（还有其他的吗？）。为了让比较更具挑战性，我想提出以下规则：

1. 输入矩阵的大小分别是2048x2048和32x32。
2. 单精度或双精度浮点数都可以。
3. 将输入矩阵转换为适当格式所花的时间不算——只算卷积这一步。
4. 用输出替换输入矩阵是可以的（有没有Python库支持这个？）
5. 直接调用常见C库的DLL是可以的——比如lapack或scalapack。
6. PyCUDA不可以。使用自定义的GPU硬件不公平。

Answer 1

我也做了一些实验。我的猜测是，SciPy的卷积计算没有使用BLAS库来加速处理。使用BLAS后，我能写出一个2D卷积，速度和MATLAB的差不多。虽然这样做需要更多的工作，但如果想要更快的速度，重新用C++编写卷积代码是个不错的选择。

下面是循环中关键的部分（请原谅我用的奇怪的()数组引用方式，这是我为MATLAB数组准备的便利类）。关键在于，你不是遍历图像，而是遍历滤波器，让BLAS来遍历图像，因为通常图像要比滤波器大得多。

for(int n = 0; n < filt.numCols; n++)
  {
    for(int m = 0; m < filt.numRows; m++)
    {
      const double filt_val = filt(filt.numRows-1-m,filt.numCols-1-n);
      for (int i =0; i < diffN; i++)
      {
        double *out_ptr = &outImage(0,i);
        const double *im_ptr = &image(m,i+n);
        cblas_daxpy(diffM,filt_val,im_ptr, 1, out_ptr,1);

      }
   }
 }

Answer 2

在我的电脑上，手动制作的圆形卷积使用快速傅里叶变换（FFT）似乎是最快的：

import numpy
x = numpy.random.random((2048, 2048)).astype(numpy.float32)
y = numpy.random.random((32, 32)).astype(numpy.float32)
z = numpy.fft.irfft2(numpy.fft.rfft2(x) * numpy.fft.rfft2(y, x.shape))

需要注意的是，这种方法可能会对靠近边缘的区域处理得和其他方法不一样，因为它是圆形卷积。

Answer 3

这真的要看你想做什么……很多时候，你并不需要一个完全通用的（也就是比较慢的）二维卷积……（也就是说，如果滤波器是可分离的，你可以用两个一维卷积来代替……这就是为什么像 scipy.ndimage.gaussian 和 scipy.ndimage.uniform 这些方法比用通用的n维卷积实现的要快得多。）

无论如何，作为一个比较：

t = timeit.timeit(stmt='ndimage.convolve(x, y, output=x)', number=1,
setup="""
import numpy as np
from scipy import ndimage
x = np.random.random((2048, 2048)).astype(np.float32)
y = np.random.random((32, 32)).astype(np.float32)
""")
print t

在我的机器上，这个过程花了6.9秒……

把它和 fftconvolve 比较一下

t = timeit.timeit(stmt="signal.fftconvolve(x, y, mode='same')", number=1,
setup="""
import numpy as np
from scipy import signal
x = np.random.random((2048, 2048)).astype(np.float32)
y = np.random.random((32, 32)).astype(np.float32)
""")
print t

这个大约花了10.8秒。不过，使用不同的输入大小时，利用傅里叶变换做卷积可能会快很多（虽然我现在想不出一个好的例子……）。