5

Matlab 中的卷积似乎比 Numpy 中的卷积快两倍。

Python 代码（在我的机器上需要 19 秒）：

import numpy as np
from scipy import ndimage
import time

img = np.ones((512,512,512))
kernel = np.ones((5,5,5))/125

start_time = time.time()
ndimage.convolve(img,kernel,mode='constant')
print "Numpy execution took ", (time.time() - start_time), "seconds"

Matlab 代码（在我的机器上需要 8.7 秒）：

img = ones(512,512,512);
kernel = ones(5,5,5) / 125;
tic
convn(img, kernel, 'same');
toc

两者给出相同的结果。

有没有办法改进 Numpy 以匹配或击败 Matlab 的性能？

有趣的是，这个因素或运行时的 ~2 差异在许多输入大小下是一致的。

4

2 回答 2

8

你在做什么具体的操作？ndimage如果您不需要一般的 Nd 卷积，可以提供许多优化。

例如，您当前的操作：

img = np.ones((512,512,512))
kernel = np.ones((5,5,5))/125
result = ndimage.convolve(img, kernel)

相当于：

img = np.ones((512,512,512))
result = ndimage.uniform_filter(img, 5)

但是，执行速度存在很大差异：

In [22]: %timeit ndimage.convolve(img, kernel)
1 loops, best of 3: 25.9 s per loop

In [23]: %timeit ndimage.uniform_filter(img, 5)
1 loops, best of 3: 8.69 s per loop

差异是由uniform_filter沿每个轴执行一维卷积而不是通用的 3D 卷积引起的。

在内核是对称的情况下，您可以进行这些简化并获得显着的加速。

我不确定convn，但如果您的输入数据符合某些标准，matlab 的函数通常会在幕后进行此类优化。Scipy 更常用的是每个函数一个算法，并希望用户知道在哪种情况下选择哪个算法。

您提到了“高斯拉普拉斯算子”过滤器。我总是对这里的术语感到困惑。

我认为你想要的ndimage功能是要么scipy.ndimage.gaussian_laplace或scipy.ndimage.gaussian_filterwith order=2（它通过高斯的二阶导数过滤）。

无论如何，两者都将操作简化为每个轴上的一维卷积，这应该会产生显着的（2-3 倍）加速。

于 2013-05-23T18:31:51.127 回答
3

不是答案；只是一个评论：

在比较性能之前，您需要确保两个函数返回相同的结果：

如果 Matlab 的 convn 返回与 Octave 的 convn 相同的结果，则 convn不同于ndimage.convolve：

octave> convn(ones(3,3), ones(2,2))
ans =

   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1

In [99]: ndimage.convolve(np.ones((3,3)), np.ones((2,2)))
Out[99]: 
array([[ 4.,  4.,  4.],
       [ 4.,  4.,  4.],
       [ 4.,  4.,  4.]])

ndimage.convolve有其他模式，'reflect','constant','nearest','mirror', 'wrap'但这些都不匹配convn的默认（“完整”）行为。

对于 2D 数组，scipy.signal.convolve2d比scipy.signal.convolve.

对于 3D 数组，scipy.signal.convolve似乎具有与以下相同的行为convn(A,B)：

octave> x = convn(ones(3,3,3), ones(2,2,2))
x =

ans(:,:,1) =

   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1

ans(:,:,2) =

   2   4   4   2
   4   8   8   4
   4   8   8   4
   2   4   4   2

ans(:,:,3) =

   2   4   4   2
   4   8   8   4
   4   8   8   4
   2   4   4   2

ans(:,:,4) =

   1   2   2   1
   2   4   4   2
   2   4   4   2
   1   2   2   1

In [109]: signal.convolve(np.ones((3,3,3), dtype='uint8'), np.ones((2,2,2), dtype='uint8'))
Out[109]: 
array([[[1, 2, 2, 1],
        [2, 4, 4, 2],
        [2, 4, 4, 2],
        [1, 2, 2, 1]],

       [[2, 4, 4, 2],
        [4, 8, 8, 4],
        [4, 8, 8, 4],
        [2, 4, 4, 2]],

       [[2, 4, 4, 2],
        [4, 8, 8, 4],
        [4, 8, 8, 4],
        [2, 4, 4, 2]],

       [[1, 2, 2, 1],
        [2, 4, 4, 2],
        [2, 4, 4, 2],
        [1, 2, 2, 1]]], dtype=uint8)

请注意，默认情况下np.ones((n,m,p))会创建一个浮点数组。Matlabones(n,m,p)似乎创建了一个整数数组。为了进行良好的比较，您应该尝试将 numpy 数组的 dtype 与 Matlab 矩阵的类型相匹配。

于 2013-05-23T17:42:25.807 回答