views:

2789

answers:

4

Is there a FFT-based 2D cross-correlation or convolution function built into scipy (or another popular library)? There are functions like these:

  • scipy.signal.correlate2d - "the direct method implemented by convolveND will be slow for large data"
  • scipy.ndimage.correlate - "The array is correlated with the given kernel using exact calculation (i.e. not FFT)."
  • scipy.fftpack.convolve.convolve, which I don't really understand, but seems wrong

Numarray had a correlate2d() function, with a 'fft=True' switch (http://structure.usc.edu/numarray/node61.html), but I guess numarray was folded into numpy, and I can't find if this function was included.

+2  A: 

I think you want the scipy.stsci package:

http://docs.scipy.org/doc/scipy/reference/stsci.html

In [30]: scipy.__version__
Out[30]: '0.7.0'

In [31]: from scipy.stsci.convolve import convolve2d, correlate2d
ars
I saw that, too, but it doesn't seem to be included in SciPy anymore?>>> import scipy.stsci.convolveTraceback (most recent call last): File "<stdin>", line 1, in <module>ImportError: No module named convolve
endolith
Hi - I pasted the output from my prompt above. What's your version?
ars
Clearly something is wrong:http://pastebin.com/mdd2bc6dGood to know it exists, though.
endolith
Werid. From your ipython prompt, I see you're using python 2.6. I have python 2.5.2. I have no idea why scipy would have a different release per version. Maybe it's easier to just re-install scipy and see if the problem persists?
ars
It works on my Windows machine with 2.6, but not on other Ubuntu machines, so it must be a packaging issue with Ubuntu. https://bugs.launchpad.net/bugs/397217
endolith
you could use correlate2d from scipy.signal: it uses more or less the same implementation technique as the stsci.convolve one (no FFT). The 2.6 problem is weird - may be related to a distutils thing.
David Cournapeau
+2  A: 

I've lost track of the status of this package in scipy, but I know we include ndimage as part of the stsci_python release package as a convenience for our users:

http://www.stsci.edu/resources/software_hardware/pyraf/stsci_python/current/download

or you should be able pull it from the repository if you prefer:

https://www.stsci.edu/svn/ssb/stsci_python/stsci_python/trunk/ndimage/

Vicki Laidler
According to SciPy docs it's not FFT-based, though, as I mentioned in the question. http://www.scipy.org/SciPyPackages/Ndimage
endolith
The convolve package is also available from the stsci_python repository. It includes the correlate2d function that has the fft=True switch that you also mentioned. https://www.stsci.edu/svn/ssb/stsci_python/stsci_python/trunk/convolve/lib/Convolve.py
Vicki Laidler
Oh! I can just import that python file directly, if I remove the reference to _correlate. The FFT correlation is all in Python. Now I've got it working. :) Thanks!
endolith
Turns out stsci is being removed from SciPy (which is why it doesn't work) and the stsci_python version is now the authoritative one, so I'm moving this to be the accepted answer.
endolith
Also, the 1-D convolve/correlate is not FFT-accelerated. :(
endolith
+2  A: 

look at scipy.signal.fftconvolve, signal.convolve and signal.correlate (there is a signal.correlate2d but it seems to return an shifted array, not centered).

magnus
I changed my accepted answer to this, as explained below http://stackoverflow.com/questions/1100100/fft-based-2d-convolution-and-correlation-in-python/1768140#1768140
endolith
+1  A: 

I found scipy.signal.fftconvolve, as also pointed out by magnus, but didn't realize at the time that it's n-dimensional. Since it's built-in and produces the right values, it seems like the ideal solution.

From Example of 2D Convolution:

In [1]: a = asarray([[ 1, 2, 3],
   ...:              [ 4, 5, 6],
   ...:              [ 7, 8, 9]])

In [2]: b = asarray([[-1,-2,-1],
   ...:              [ 0, 0, 0],
   ...:              [ 1, 2, 1]])

In [3]: scipy.signal.fftconvolve(a, b, mode = 'same')
Out[3]: 
array([[-13., -20., -17.],
       [-18., -24., -18.],
       [ 13.,  20.,  17.]])

Correct! The STSCI version, on the other hand, requires some extra work to make the boundaries correct?

In [4]: stsci.convolve2d(a, b, fft = True)
Out[4]: 
array([[-12., -12., -12.],
       [-24., -24., -24.],
       [-12., -12., -12.]])

The STSCI method also requires compiling, which I was unsuccessful with (I just commented out the non-python parts), has some bugs like this and modifying the inputs ([1, 2] becomes [[1, 2]]), etc. So I changed my accepted answer to the built-in fftconvolve() function. Correlation, of course, is the same thing as convolution, but with one input reversed:

In [5]: a
Out[5]: 
array([[3, 0, 0],
       [2, 0, 0],
       [1, 0, 0]])

In [6]: b
Out[6]: 
array([[3, 2, 1],
       [0, 0, 0],
       [0, 0, 0]])

In [7]: scipy.signal.fftconvolve(a, b[::-1, ::-1])
Out[7]: 
array([[ 0., -0.,  0.,  0.,  0.],
       [ 0., -0.,  0.,  0.,  0.],
       [ 3.,  6.,  9.,  0.,  0.],
       [ 2.,  4.,  6.,  0.,  0.],
       [ 1.,  2.,  3.,  0.,  0.]])

In [8]: scipy.signal.correlate2d(a, b)
Out[8]: 
array([[0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [3, 6, 9, 0, 0],
       [2, 4, 6, 0, 0],
       [1, 2, 3, 0, 0]])

and the latest revision has been sped up by using power-of-two sizes internally.

endolith