24

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 an fft=True switch, but I guess numarray was folded into numpy, and I can't find if this function was included.

ali_m
  • 71,714
  • 23
  • 223
  • 298
endolith
  • 25,479
  • 34
  • 128
  • 192
  • 1
    note that using exact calculation (no FFT) is exactly the same as saying it is slow :) More exactly, the FFT-based method will be much faster if you have a signal and a kernel of approximately the same size (if the kernel is much smaller than the input, then FFT may actually be slower than the direct computation). – David Cournapeau Jul 26 '09 at 08:18
  • Ideally, the FFT algorithm would automatically take care of zero-padding things to the right size for best speed. – endolith Aug 17 '09 at 18:26
  • 1
    Oh you're not talking about zero padding, you're talking about matching a 5x5 image with a 2000x2000 image. Why can't the algorithm just guess whether the FFT would be more efficient and do it whichever way is faster? – endolith Aug 20 '09 at 23:04
  • 2
    scipy has an fftconvolve function http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html#scipy.signal.fftconvolve – endolith Sep 01 '09 at 17:48

6 Answers6

25

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 (and then I sped it up more by using real FFT for real input and using 5-smooth lengths instead of powers of 2 :D ).

Community
  • 1
  • 1
endolith
  • 25,479
  • 34
  • 128
  • 192
8

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).

  • 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 Nov 20 '09 at 15:16
  • Hi, just an unrelated query. How long does correlate2d take to run for 2-d arrays of shape 2048x2048? – lordparthurnaax May 16 '19 at 20:34
5

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
  • 120,335
  • 23
  • 147
  • 134
  • I saw that, too, but it doesn't seem to be included in SciPy anymore? >>> import scipy.stsci.convolve Traceback (most recent call last): File "", line 1, in ImportError: No module named convolve – endolith Jul 08 '09 at 21:01
  • Hi - I pasted the output from my prompt above. What's your version? – ars Jul 08 '09 at 21:08
  • Clearly something is wrong: http://pastebin.com/mdd2bc6d Good to know it exists, though. – endolith Jul 08 '09 at 23:16
  • 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 Jul 08 '09 at 23:27
  • 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 Jul 13 '09 at 18:53
  • 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 Jul 26 '09 at 08:19
4

I wrote a cross-correlation/convolution wrapper that takes care of padding & nans and includes a simple smooth wrapper here. It's not a popular package, but it also has no dependencies besides numpy (or fftw for faster ffts).

I've also implemented an FFT speed testing code here in case anyone's interested. It shows - surprisingly - that numpy's fft is faster than scipy's, at least on my machine.

EDIT: moved code to N-dimensional version here

keflavich
  • 18,278
  • 20
  • 86
  • 118
3

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
  • 3,415
  • 1
  • 20
  • 17
  • According to SciPy docs it's not FFT-based, though, as I mentioned in the question. http://www.scipy.org/SciPyPackages/Ndimage – endolith Aug 15 '09 at 19:13
  • 1
    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 Aug 16 '09 at 05:14
  • 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 Aug 17 '09 at 18:28
  • 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 Aug 21 '09 at 00:33
  • 1
    Also, the 1-D convolve/correlate is not FFT-accelerated. :( – endolith Sep 01 '09 at 17:44
2

Note that there is scipy.signal.oaconvolve since Scipy 1.4, which uses the Overlap-add method.

L. F. Ant
  • 21
  • 3