Python 中最快的 2D 卷积或图像过滤器答案

【问题标题】：Fastest 2D convolution or image filter in PythonPython 中最快的 2D 卷积或图像过滤器
【发布时间】：2011-08-08 07:48:22
【问题描述】：

一些用户询问过 numpy 或 scipy 中图像卷积的速度或内存消耗 [1、2、3、4]。从回复和我使用 Numpy 的经验来看，我认为这可能是 numpy 与 Matlab 或 IDL 相比的主要缺点。

到目前为止，没有一个答案解决了整个问题，所以这里是：“在 Python 中计算 2D 卷积的最快方法是什么？”常见的 python 模块是公平的游戏：numpy、scipy 和 PIL（其他？）。为了具有挑战性的比较，我想提出以下规则：

输入矩阵分别为 2048x2048 和 32x32。
单精度或双精度浮点均可接受。
不计算将输入矩阵转换为适当格式所花费的时间 - 仅计算卷积步骤。
用您的输出替换输入矩阵是可以接受的（有任何 python 库支持吗？）
对常用 C 库的直接 DLL 调用没问题 -- lapack 或 scalapack
PyCUDA 是正确的。使用自定义 GPU 硬件是不公平的。

【问题讨论】：

“用你的输出替换输入矩阵是可以接受的（任何 python 库都支持吗？）”对于它的价值，大多数 numpy 和 scipy 函数都可以......
我在 convolve 的文档中没有看到任何提及：docs.scipy.org/doc/numpy/reference/generated/… 我错过了什么吗？
numpy 的 convolve 不支持，scipy.ndimage.convolve 支持。 scipy.org/SciPyPackages/Ndimage 此外，大多数 numpy 函数（例如 sqrt、mul、add）都带有一个 out 参数。您可以使用np.sqrt(x, x) 将 sqrt 就地执行。

标签： python optimization numpy python-imaging-library scipy

【解决方案1】：

这真的取决于你想做什么......很多时候，你不需要一个完全通用的（阅读：更慢的）2D卷积......（即如果过滤器是可分离的，你使用两个1D 卷积代替...这就是为什么各种 scipy.ndimage.gaussian、scipy.ndimage.uniform 比作为通用 nD 卷积实现的相同事物快得多的原因。）

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

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 进行卷积可能会快得多（尽管目前我似乎无法提出一个很好的例子......）。

【讨论】：

谢谢乔。这是对我一直使用的 convolve 函数的重大改进（我认为它只是 numpy.convolve）。它消耗了大量的内存并且运行缓慢（可能是结果）。我希望得到更多的参与，但也许我太乐观了。
对于那些感兴趣的人。我做了这个比较（OS X 10.10 Macbook Air）比原来的帖子晚了 5 年。 signal.fftconvolve 大约需要 .9 秒！ ndimage.convolve 大约需要 8 秒。 signal.fftconvolve 的底层显然已经进行了巨大的改进。

【解决方案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))

请注意，这可能会以不同于其他方式的方式处理靠近边缘的区域，因为它是循环卷积。

【讨论】：

无限快
您如何将相同的概念应用于“相同”的卷积？

【解决方案3】：

我也对此做了一些实验。我的猜测是 SciPy 卷积不使用 BLAS 库来加速计算。使用 BLAS，我能够编写速度与 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);

      }
   }
 }

【讨论】：

【解决方案4】：

我一直在尝试提高应用程序中的卷积速度，并且我一直在使用 signal.correlate，它恰好比 signal.correlate2d 慢大约 20 倍，我的输入矩阵更小（27x27 and 5x5）。截至 2018 年，这是我在我的机器（Dell Inspiron 13，Core i5）上观察到的实际问题中指定矩阵的结果。

OpenCV 做得最好，但需要注意的是它没有提供“模式”选项。输入和输出的大小相同。

>>> img= np.random.rand(2048,2048)
>>> kernel = np.ones((32,32), dtype=np.float)
>>> t1= time.time();dst1 = cv2.filter2D(img,-1,kernel);print(time.time()-t1)
0.208490133286
>>> t1= time.time();dst2 = signal.correlate(img,kernel,mode='valid',method='fft');print(time.time()-t1)
0.582989931107
>>> t1= time.time();dst3 = signal.convolve2d(img,kernel,mode='valid');print(time.time()-t1)
11.2672450542
>>> t1= time.time();dst4 = signal.correlate2d(img,kernel,mode='valid');print(time.time()-t1)
11.2443971634
>>> t1= time.time();dst5 = signal.fftconvolve(img,kernel,mode='valid');print(time.time()-t1)
0.581533193588

【讨论】：

【解决方案5】：

Scipy 有一个函数fftconvolve，可用于一维和二维信号。

from scipy import signal
from scipy import misc
import numpy as np
import matplotlib.pyplot as plt

face = misc.face(gray=True)
kernel = np.outer(signal.gaussian(70, 8), signal.gaussian(70, 8))
blurred = signal.fftconvolve(face, kernel, mode='same')

fig, (ax_orig, ax_kernel, ax_blurred) = plt.subplots(3, 1, figsize=(6, 15))
ax_orig.imshow(face, cmap='gray')
ax_orig.set_title('Original')
ax_orig.set_axis_off()
ax_kernel.imshow(kernel, cmap='gray')
ax_kernel.set_title('Gaussian kernel')
ax_kernel.set_axis_off()
ax_blurred.imshow(blurred, cmap='gray')
ax_blurred.set_title('Blurred')
ax_blurred.set_axis_off()
fig.show()

【讨论】：