Computing cross-correlation function?

Question

In R, I am using ccf or acf to compute the pair-wise cross-correlation function so that I can find out which shift gives me the maximum value. From the looks of it, R gives me a normalized sequence of values. Is there something similar in Python's scipy or am I supposed to do it using the fft module? Currently, I am doing it as follows:

xcorr = lambda x,y : irfft(rfft(x)*rfft(y[::-1]))
x = numpy.array([0,0,1,1])
y = numpy.array([1,1,0,0])
print xcorr(x,y)

Answer 1

To cross-correlate 1d arrays use numpy.correlate.

For 2d arrays, use scipy.signal.correlate2d.

There is also scipy.stsci.convolve.correlate2d.

There is also matplotlib.pyplot.xcorr which is based on numpy.correlate.

See this post on the SciPy mailing list for some links to different implementations.

Edit: @user333700 added a link to the SciPy ticket for this issue in a comment.

Answer 2

If you are looking for a rapid, normalized cross correlation in either one or two dimensions I would recommend the openCV library (see http://opencv.willowgarage.com/wiki/ http://opencv.org/). The cross-correlation code maintained by this group is the fastest you will find, and it will be normalized (results between -1 and 1).

While this is a C++ library the code is maintained with CMake and has python bindings so that access to the cross correlation functions is convenient. OpenCV also plays nicely with numpy. If I wanted to compute a 2-D cross-correlation starting from numpy arrays I could do it as follows.

import numpy
import cv

#Create a random template and place it in a larger image
templateNp = numpy.random.random( (100,100) )
image = numpy.random.random( (400,400) )
image[:100, :100] = templateNp

#create a numpy array for storing result
resultNp = numpy.zeros( (301, 301) )

#convert from numpy format to openCV format
templateCv = cv.fromarray(numpy.float32(template))
imageCv = cv.fromarray(numpy.float32(image))
resultCv =  cv.fromarray(numpy.float32(resultNp))

#perform cross correlation
cv.MatchTemplate(templateCv, imageCv, resultCv, cv.CV_TM_CCORR_NORMED)

#convert result back to numpy array
resultNp = np.asarray(resultCv)

For just a 1-D cross-correlation create a 2-D array with shape equal to (N, 1 ). Though there is some extra code involved to convert to an openCV format the speed-up over scipy is quite impressive.

Answer 3

I just finished writing my own optimised implementation of normalized cross-correlation for N-dimensional arrays. You can get it from here.

It will calculate cross-correlation either directly, using scipy.ndimage.correlate, or in the frequency domain, using scipy.fftpack.fftn/ifftn depending on whichever will be quickest.

Answer 4

For 1D array, numpy.correlate is faster than scipy.signal.correlate, under different sizes, I see a consistent 5x peformance gain using numpy.correlate. When two arrays are of similar size (the bright line connecting the diagonal), the performance difference is even more outstanding (50x +).

# a simple benchmark
res = []
for x in range(1, 1000):
    list_x = []
    for y in range(1, 1000): 

        # generate different sizes of series to compare
        l1 = np.random.choice(range(1, 100), size=x)
        l2 = np.random.choice(range(1, 100), size=y)

        time_start = datetime.now()
        np.correlate(a=l1, v=l2)
        t_np = datetime.now() - time_start

        time_start = datetime.now()
        scipy.signal.correlate(in1=l1, in2=l2)
        t_scipy = datetime.now() - time_start

        list_x.append(t_scipy / t_np)
    res.append(list_x)
plt.imshow(np.matrix(res))

As default, scipy.signal.correlate calculates a few extra numbers by padding and that might explained the performance difference.

>> l1 = [1,2,3,2,1,2,3]
>> l2 = [1,2,3]
>> print(numpy.correlate(a=l1, v=l2))
>> print(scipy.signal.correlate(in1=l1, in2=l2))

[14 14 10 10 14]
[ 3  8 14 14 10 10 14  8  3]  # the first 3 is [0,0,1]dot[1,2,3]