13

As far as I have seen, these methods are both implemented as C functions in the respective DLLs, and it appears that the ndimage version is faster (neither implementation uses parallelized code, like calls to blas or MKL).

Also, when I tried to check that they return the same results by running the following code, the assertion of equality failed. I couldn't figure out from the documentation what exactly the functional differences between the two methods should be (the documentation isn't very clear about what 0 means relative to the location of the kernel's origin either; from examples, I deduced it's in the center, but I might be wrong).

from numpy import random, allclose
from scipy.ndimage.filters import convolve as convolveim
from scipy.signal import convolve as convolvesig

a = random.random((100, 100, 100))
b = random.random((10,10,10))

conv1 = convolveim(a,b, mode = 'constant')
conv2 = convolvesig(a,b, mode = 'same')

assert(allclose(conv1,conv2))

Thanks!

Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
bbudescu
  • 313
  • 1
  • 3
  • 9
  • 1
    I do not know the implementations, but probably the implementation from ndimage uses the Convolution Theorem, i.e., convolution is equal to multiplication in Fourier space. This is what scipy.signal.fftconvolve does. But also when using this method instead of convolve, the assertion fails. – Thorsten Kranz Jun 03 '13 at 11:39
  • 1
    It seems to me that both of them are using direct implementations. I beleieve that the code that gets called in the end is [here](https://github.com/scipy/scipy/blob/master/scipy/signal/correlate_nd.c.src#L105) and [here](https://github.com/scipy/scipy/blob/master/scipy/ndimage/src/ni_filters.c#L132), respectively. I tried calling the methods on integer arrays to rule out rounding errors etc., but the assertion fails as well. – bbudescu Jun 03 '13 at 11:57
  • 1
    The assertion fails on 2d arrays as well, but only if the kernel is greater than a certain size e.g. `a=random.random_integers(0,10,(100, 100));b=random.random_integers(0,10,(7, 7))` doesn't fail, but when `b=random.random_integers(0,10,(8, 8))`, it does. Any thoughts? – bbudescu Jun 03 '13 at 12:01
  • [Great answer here](http://stackoverflow.com/a/16121975/3079302). – iled Feb 24 '17 at 04:48

2 Answers2

13

The two functions have different conventions for dealing with the boundary. To make your calls functionally the same, add the argument origin=-1 or origin=(-1,-1,-1) to the call to convolveim:

In [46]: a = random.random((100,100,100))

In [47]: b = random.random((10,10,10))

In [48]: c1 = convolveim(a, b, mode='constant', origin=-1)

In [49]: c2 = convolvesig(a, b, mode='same')

In [50]: allclose(c1,c2)
Out[50]: True

Only shift the origin when the dimensions of b are even. When they are odd, the functions agree when the default origin=0 is used:

In [88]: b = random.random((11,11,11))

In [89]: c1 = convolveim(a, b, mode='constant')

In [90]: c2 = convolvesig(a, b, mode='same')

In [91]: allclose(c1,c2)
Out[91]: True
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • If a and b have different sizes (modulo 2), e.g. if a has size (100,100) and b has size (11,11), the above solution is incorrect. In this case the origin should be (0,0) for the two functions to agree – Jens Munk Jul 03 '14 at 08:35
  • @JensMunk What exactly did you test? If I convolve a (1000,1000) array with a (11,11) array with `origin=0`, the results are very close (essentially the same exact some numerical precision). This is what the answer says to do if `b` is odd. – user3731622 Mar 09 '17 at 20:43
  • That's right. If they have the same size modulo 2 the origin should be (1,1) as far as I remember. – Jens Munk Mar 10 '17 at 03:41
5

There is a very important difference. The implementation in the image package seems to be the typical restricted version used in image processing to achieve the 'same' size of the image after the convolution. So it coincides with the 'same' option in the signal processing package if we use the mode='constant', as in the examples above. The signal processing package seems to implement the real strict definition of the convolution operator. Perhaps for this reason it is slower. Find enclosed some examples with completely different results.

In [13]: a=array([[1,2,1]])
In [14]: b=array([[1],[2],[1]])

In [17]: convolveim(a,b)
Out[17]: array([[4, 8, 4]])

In [18]: convolveim(b,a)
Out[18]: 
array([[4],
       [8],
       [4]])

In [19]: convolvesig(a,b)
Out[19]: 
array([[1, 2, 1],
       [2, 4, 2],
       [1, 2, 1]])

In [20]: convolvesig(b,a)
Out[20]: 
array([[1, 2, 1],
       [2, 4, 2],
       [1, 2, 1]])

Note that the signal processing package implementation is conmutative, as expected for a correct convolution. However, the implementation in the image package is not and it provides a solution with the 'same' dimensions as the first parameter.

kapitan_tan
  • 51
  • 1
  • 2