Phase Correlation in Opencv

Phase Correlation in OpenCV

Phase Correlation is a method to check the similarity of two images with equal size. It can be used for template matching, object tracking, motion estimation, etc

In this article, I explain the basics and the code to perform Phase Correlation in OpenCV, divided into several chapters:

1. The Basics
Intoduction to Phase Correlation.

2. Main Code
Line-by-line explanation of main.c.

3. Phase Correlation Function
Line-by-line explanation of phase_correlation.c.

4. Results
The results from some experiment with the code.

The source code used in this tutorial is available for you to download. The package contains all necessary files and some installation notes.

2. The Basics

To obtain the Phase Correlation of two images, perform these steps:

1. Load two images, f and g.

2. Perform FFT on each image, resulting in F and G

3. Obtain the cross power spectrum using this formula:

where is the complex conjugate of G.

4. Obtain the phase correlation by performing IFFT on R

The result is a 2D array with each element has a value between 0 to 1. The location of the highest value corresponds with the object translation movement from image 1 to image 2.

3. Main Code

The main code performs a very simple task. It loads the reference and template images, call phase_correlation() function and display the result.

Note: In the code listings below I removed all error checking lines to make things simple. Remember that you should always have some error checking in your programs.

First, load the template and the reference image. Since Phase Correlation computes 2 equal images, the width and the height of both images need to be checked.

Listing 1: Load reference and template image

/* load reference and template image */

2. IplImage *ref = cvLoadImage

IplImage *tpl = cvLoadImage( argv[2], CV_LOAD_IMAGE_GRAYSCALE );

5. /* both images' size should be equal */

if( ( tpl->width != ref->width ) || ( tpl->height != ref->height ) ) {
fprintf( stderr, "Both images must have equal width and height!\n" );
return 1;
}

10.

Next create a new image, poc, for the Phase Correlation result. Then call phase_correlation() function with tpl, ref and poc passed as parameters.

Listing 2: Perform Phase Correlation

11. /* create a new image, to store phase correlation result */

12. IplImage *poc = cvCreateImage( cvSize( tpl->width, tpl->height ),

13. IPL_DEPTH_64F, 1 );

14.

15. /* get phase correlation of input images */

16. phase_correlation( ref, tpl, poc );

17.

Now that poc contains the result from Phase Correlation, find its maximum value and the location using cvMinMaxLoc. The result then displayed in the command line.

Listing 3: Get maximum value

18. /* find the maximum value and its location */

19. CvPoint minloc, maxloc;

20. double minval, maxval;

21. cvMinMaxLoc( poc, &minval, &maxval, &minloc, &maxloc, 0 );

22.

23. /* print it */

24. fprintf( stdout,

25. "Maxval at (%d, %d) = %2.4f\n", maxloc.x, maxloc.y, maxval );

The last step, display the images and free memory.

Listing 4: Display images and free memory

27. /* display images and free memory */

28. cvNamedWindow( "tpl", CV_WINDOW_AUTOSIZE );

29. cvNamedWindow( "ref", CV_WINDOW_AUTOSIZE );

30. cvShowImage( "tpl", tpl );

31. cvShowImage( "ref", ref );

32.

33. cvWaitKey( 0 );

34.

35. cvDestroyWindow( "tpl" );

36. cvDestroyWindow( "ref" );

37. cvReleaseImage( &tpl );

38. cvReleaseImage( &ref );

39. cvReleaseImage( &poc );

In the next chapter I will explain the phase_correlation() function, line by line

4. Phase Correlation Function

This function performs Phase Correlation operation. It takes 3 parameters:

1. IplImage *ref: The reference image

2. IplImage *tpl: The template image

3. IplImage *poc: Empty image to store the result

To obtain the DFT of the images, I use the FFTW library rather than OpenCV's DFT functions. This will increase the speed since FFTW computes FFT very fast.

First get the width and the height of the image. These properties are needed for the rest of the function. Also setup the pointers to copy the images to FFTW input and vice versa.

Listing 5: Phase Correlation Function

1. void phase_correlation( IplImage *ref, IplImage *tpl, IplImage *poc )

2. {

3. int i, j, k;

4. double tmp;

6. /* get image properties */

7. int width = ref->width;

8. int height = ref->height;

9. int step = ref->widthStep;

10. int fft_size = width * height;

11.

12. /* setup pointers to images */

13. uchar *ref_data = ( uchar* ) ref->imageData;

14. uchar *tpl_data = ( uchar* ) tpl->imageData;

15. double *poc_data = ( double* )poc->imageData;

16.

Note that the function above assumes that the three images have equal width and height. Next initialize some FFTW input/output arrays and plans, load the reference image to img1 and the template image to img2.

Listing 6: Initialize FFTW input and output arrays

17. /* allocate FFTW input and output arrays */

18. fftw_complex *img1 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

19. fftw_complex *img2 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

20. fftw_complex *res = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

21.

22. /* setup FFTW plans */

23. fftw_plan fft_img1 = fftw_plan_dft_1d( width * height, img1, img1, FFTW_FORWARD, FFTW_ESTIMATE );

24. fftw_plan fft_img2 = fftw_plan_dft_1d( width * height, img2, img2, FFTW_FORWARD, FFTW_ESTIMATE );

25. fftw_plan ifft_res = fftw_plan_dft_1d( width * height, res, res, FFTW_BACKWARD, FFTW_ESTIMATE );

26.

27. /* load images' data to FFTW input */

28. for( i = 0, k = 0 ; i < height ; i++ ) {

29. for( j = 0 ; j < width ; j++, k++ ) {

30. img1[k][0] = ( double )ref_data[i * step + j];

31. img1[k][1] = 0.0;

32.

33. img2[k][0] = ( double )tpl_data[i * step + j];

34. img2[k][1] = 0.0;

35. }

36. }

37.

In the listing above res will hold the end result of the phase correlation operation. Next is the heart of this function. It performs Phase Correlation and store the result to poc.

Listing 7: Compute phase correlation

38. /* obtain the FFT of img1 */

39. fftw_execute( fft_img1 );

40.

41. /* obtain the FFT of img2 */

42. fftw_execute( fft_img2 );

43.

44. /* obtain the cross power spectrum */

45. for( i = 0; i < fft_size ; i++ ) {

46. res[i][0] = ( img2[i][0] * img1[i][0] ) - ( img2[i][1] * ( -img1[i][1] ) );

47. res[i][1] = ( img2[i][0] * ( -img1[i][1] ) ) + ( img2[i][1] * img1[i][0] );

48.

49. tmp = sqrt( pow( res[i][0], 2.0 ) + pow( res[i][1], 2.0 ) );

50.

51. res[i][0] /= tmp;

52. res[i][1] /= tmp;

53. }

54.

55. /* obtain the phase correlation array */

56. fftw_execute(ifft_res);

57.

58. /* normalize and copy to result image */

59. for( i = 0 ; i < fft_size ; i++ ) {

60. poc_data[i] = res[i][0] / ( double )fft_size;

61. }

62.

Note that we have to normalize the IFFT result since FFTW computes an unnormalized DFT. Now all we have to do is just free memory and return to main function.

Listing 8: Free Memory

63. /* deallocate FFTW arrays and plans */

64. fftw_destroy_plan( fft_img1 );

65. fftw_destroy_plan( fft_img2 );

66. fftw_destroy_plan( ifft_res );

67. fftw_free( img1 );

68. fftw_free( img2 );

69. fftw_free( res );

70. }

The Phase Correlation result is stored in poc. The main code should find the maximum value and its location from this result.

1. IplImage *ref: The reference image

2. IplImage *tpl: The template image

3. IplImage *poc: Empty image to store the result

To obtain the DFT of the images, I use the FFTW library rather than OpenCV's DFT functions. This will increase the speed since FFTW computes FFT very fast.

First get the width and the height of the image. These properties are needed for the rest of the function. Also setup the pointers to copy the images to FFTW input and vice versa.

Listing 5: Phase Correlation Function

1. void phase_correlation( IplImage *ref, IplImage *tpl, IplImage *poc )

2. {

3. int i, j, k;

4. double tmp;

6. /* get image properties */

7. int width = ref->width;

8. int height = ref->height;

9. int step = ref->widthStep;

10. int fft_size = width * height;

11.

12. /* setup pointers to images */

13. uchar *ref_data = ( uchar* ) ref->imageData;

14. uchar *tpl_data = ( uchar* ) tpl->imageData;

15. double *poc_data = ( double* )poc->imageData;

16.

Listing 6: Initialize FFTW input and output arrays

17. /* allocate FFTW input and output arrays */

18. fftw_complex *img1 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

19. fftw_complex *img2 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

20. fftw_complex *res = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );

21.

22. /* setup FFTW plans */

23. fftw_plan fft_img1 = fftw_plan_dft_1d( width * height, img1, img1, FFTW_FORWARD, FFTW_ESTIMATE );

24. fftw_plan fft_img2 = fftw_plan_dft_1d( width * height, img2, img2, FFTW_FORWARD, FFTW_ESTIMATE );

25. fftw_plan ifft_res = fftw_plan_dft_1d( width * height, res, res, FFTW_BACKWARD, FFTW_ESTIMATE );

26.

27. /* load images' data to FFTW input */

28. for( i = 0, k = 0 ; i < height ; i++ ) {

29. for( j = 0 ; j < width ; j++, k++ ) {

30. img1[k][0] = ( double )ref_data[i * step + j];

31. img1[k][1] = 0.0;

32.

33. img2[k][0] = ( double )tpl_data[i * step + j];

34. img2[k][1] = 0.0;

35. }

36. }

37.

In the listing above res will hold the end result of the phase correlation operation. Next is the heart of this function. It performs Phase Correlation and store the result to poc.

Listing 7: Compute phase correlation

38. /* obtain the FFT of img1 */

39. fftw_execute( fft_img1 );

40.

41. /* obtain the FFT of img2 */

42. fftw_execute( fft_img2 );

43.

44. /* obtain the cross power spectrum */

45. for( i = 0; i < fft_size ; i++ ) {

46. res[i][0] = ( img2[i][0] * img1[i][0] ) - ( img2[i][1] * ( -img1[i][1] ) );

47. res[i][1] = ( img2[i][0] * ( -img1[i][1] ) ) + ( img2[i][1] * img1[i][0] );

48.

49. tmp = sqrt( pow( res[i][0], 2.0 ) + pow( res[i][1], 2.0 ) );

50.

51. res[i][0] /= tmp;

52. res[i][1] /= tmp;

53. }

54.

55. /* obtain the phase correlation array */

56. fftw_execute(ifft_res);

57.

58. /* normalize and copy to result image */

59. for( i = 0 ; i < fft_size ; i++ ) {

60. poc_data[i] = res[i][0] / ( double )fft_size;

61. }

62.

Note that we have to normalize the IFFT result since FFTW computes an unnormalized DFT. Now all we have to do is just free memory and return to main function.

Listing 8: Free Memory

63. /* deallocate FFTW arrays and plans */

64. fftw_destroy_plan( fft_img1 );

65. fftw_destroy_plan( fft_img2 );

66. fftw_destroy_plan( ifft_res );

67. fftw_free( img1 );

68. fftw_free( img2 );

69. fftw_free( res );

70. }

The Phase Correlation result is stored in poc. The main code should find the maximum value and its location from this result.

5. Results

Here are some results I've got when computing phase correlation using this code.

The first result shows the maximum value is 1.0 at (0, 0). This means both images are identical.

The second result shows the maximum value is 0.2729 at (0, 0). This means there is no translation from image 1 to image 2, but image 2 has some noise in it.

The last result shows the maximum value is 0.6344 at (6, 15). This means image 1 has a translative movement 6 pixels to the right and 15 pixels to the bottom.

Table 1. Phase Correlation Results

# / Image 1 / Image 2 / Result /
1. / /
(Image 2 is identical
with image 1) / Max. value = 1.0
in (0,0)
2. / /
(image 1 with some noise) / Max. value = 0.2729
in (0,0)
3. / /
(object moved 6 pixels to
the right and 15 pixels to
the bottom) / Max. value = 0.6344
in (6,15)

6. Conclusion

This article covered Phase Correlation function in OpenCV. The code makes use FFTW library to obtain the DFTs rather than using OpenCV's DFT functions. This will increase the speed since FFTW computes FFT very fast.