Python中最快的二维卷积或图像滤波

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. 皮库达马上就出来了。使用自定义GPU硬件是不公平的。

Answer 1

在我的机器上，使用FFTs的手工圆形卷积似乎很快：

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 2

我也做了一些实验。我的猜测是SciPy卷积没有使用BLAS库来加速计算。使用BLAS，我能够编码一个与Matlab的速度相当的2D卷积，这是更多的工作，但你最好的办法是在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 3

这真的取决于你想做什么。。。很多时候，你不需要完全通用的二维卷积。。。（即，如果滤波器是可分离的，则使用两个一维卷积代替。。。这就是为什么不同的scipy.ndimage.gaussian，scipy.ndimage.uniform比作为普通n-D卷积实现的同一事物要快得多。）

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

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秒。但是，对于不同的输入大小，使用fft进行卷积可以更快一些（尽管目前我似乎还不能给出一个很好的例子…）。