最后是代码：
class Peak
{
public:
CvPoint pt;
double maxval;
};
Peak old_opencv_FFT(IplImage* src,IplImage* temp)
{
CvSize imgSize = cvSize(src->width, src->height);
// Allocate floating point frames used for DFT (real, imaginary, and complex)
IplImage* realInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 1 );
IplImage* imaginaryInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 1 );
IplImage* complexInput = cvCreateImage( imgSize, IPL_DEPTH_64F, 2 );
int nDFTHeight= cvGetOptimalDFTSize( imgSize.height );
int nDFTWidth= cvGetOptimalDFTSize( imgSize.width );
CvMat* src_DFT = cvCreateMat( nDFTHeight, nDFTWidth, CV_64FC2 );
CvMat* temp_DFT = cvCreateMat( nDFTHeight, nDFTWidth, CV_64FC2 );
CvSize dftSize = cvSize(nDFTWidth, nDFTHeight);
IplImage* imageRe = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
IplImage* imageIm = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
IplImage* imageImMag = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
IplImage* imageMag = cvCreateImage( dftSize, IPL_DEPTH_64F, 1 );
CvMat tmp;
// Processing of src
cvScale(src,realInput,1.0,0);
cvZero(imaginaryInput);
cvMerge(realInput,imaginaryInput,NULL,NULL,complexInput);
cvGetSubRect(src_DFT,&tmp,cvRect(0,0,src->width,src->height));
cvCopy(complexInput,&tmp,NULL);
if (src_DFT->cols>src->width)
{
cvGetSubRect(src_DFT,&tmp,cvRect(src->width,0,src_DFT->cols-src->width,src->height));
cvZero(&tmp);
}
cvDFT(src_DFT,src_DFT,CV_DXT_FORWARD,complexInput->height);
cvSplit(src_DFT,imageRe,imageIm,0,0);
// Processing of temp
cvScale(temp,realInput,1.0,0);
cvMerge(realInput,imaginaryInput,NULL,NULL,complexInput);
cvGetSubRect(temp_DFT,&tmp,cvRect(0,0,temp->width,temp->height));
cvCopy(complexInput,&tmp,NULL);
if (temp_DFT->cols>temp->width)
{
cvGetSubRect(temp_DFT,&tmp,cvRect(temp->width,0,temp_DFT->cols-temp->width,temp->height));
cvZero( &tmp );
}
cvDFT(temp_DFT,temp_DFT,CV_DXT_FORWARD,complexInput->height);
// Multiply spectrums of the scene and the model (use CV_DXT_MUL_CONJ to get correlation instead of convolution)
cvMulSpectrums(src_DFT,temp_DFT,src_DFT,CV_DXT_MUL_CONJ);
// Split Fourier in real and imaginary parts
cvSplit(src_DFT,imageRe,imageIm,0,0);
// Compute the magnitude of the spectrum components: Mag = sqrt(Re^2 + Im^2)
cvPow( imageRe, imageMag, 2.0 );
cvPow( imageIm, imageImMag, 2.0 );
cvAdd( imageMag, imageImMag, imageMag, NULL );
cvPow( imageMag, imageMag, 0.5 );
// Normalize correlation (Divide real and imaginary components by magnitude)
cvDiv(imageRe,imageMag,imageRe,1.0);
cvDiv(imageIm,imageMag,imageIm,1.0);
cvMerge(imageRe,imageIm,NULL,NULL,src_DFT);
// inverse dft
cvDFT( src_DFT, src_DFT, CV_DXT_INVERSE_SCALE, complexInput->height );
cvSplit( src_DFT, imageRe, imageIm, 0, 0 );
double minval = 0.0;
double maxval = 0.0;
CvPoint minloc;
CvPoint maxloc;
cvMinMaxLoc(imageRe,&minval,&maxval,&minloc,&maxloc,NULL);
int x=maxloc.x; // log range
//if (x>(imageRe->width/2))
// x = x-imageRe->width; // positive or negative values
int y=maxloc.y; // angle
//if (y>(imageRe->height/2))
// y = y-imageRe->height; // positive or negative values
Peak pk;
pk.maxval= maxval;
pk.pt=cvPoint(x,y);
return pk;
}
void phase_correlation2D( IplImage* src, IplImage *tpl, IplImage *poc )
{
int i, j, k;
double tmp;
/* get image properties */
int width = src->width;
int height = src->height;
int step = src->widthStep;
int fft_size = width * height;
/* setup pointers to images */
uchar *src_data = ( uchar* ) src->imageData;
uchar *tpl_data = ( uchar* ) tpl->imageData;
double *poc_data = ( double* )poc->imageData;
/* allocate FFTW input and output arrays */
fftw_complex *img1 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
fftw_complex *img2 = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
fftw_complex *res = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
/* setup FFTW plans */
fftw_plan fft_img1 = fftw_plan_dft_2d( height ,width, img1, img1, FFTW_FORWARD, FFTW_ESTIMATE );
fftw_plan fft_img2 = fftw_plan_dft_2d( height ,width, img2, img2, FFTW_FORWARD, FFTW_ESTIMATE );
fftw_plan ifft_res = fftw_plan_dft_2d( height ,width, res, res, FFTW_BACKWARD, FFTW_ESTIMATE );
/* load images' data to FFTW input */
for( i = 0, k = 0 ; i < height ; i++ ) {
for( j = 0 ; j < width ; j++, k++ ) {
img1[k][0] = ( double )src_data[i * step + j];
img1[k][1] = 0.0;
img2[k][0] = ( double )tpl_data[i * step + j];
img2[k][1] = 0.0;
}
}
///* Hamming window */
//double omega = 2.0*M_PI/(fft_size-1);
//double A= 0.54;
//double B= 0.46;
//for(i=0,k=0;i<height;i++)
//{
// for(j=0;j<width;j++,k++)
// {
// img1[k][0]= (img1[k][0])*(A-B*cos(omega*k));
// img2[k][0]= (img2[k][0])*(A-B*cos(omega*k));
// }
//}
/* obtain the FFT of img1 */
fftw_execute( fft_img1 );
/* obtain the FFT of img2 */
fftw_execute( fft_img2 );
/* obtain the cross power spectrum */
for( i = 0; i < fft_size ; i++ ) {
res[i][0] = ( img2[i][0] * img1[i][0] ) - ( img2[i][1] * ( -img1[i][1] ) );
res[i][1] = ( img2[i][0] * ( -img1[i][1] ) ) + ( img2[i][1] * img1[i][0] );
tmp = sqrt( pow( res[i][0], 2.0 ) + pow( res[i][1], 2.0 ) );
res[i][0] /= tmp;
res[i][1] /= tmp;
}
/* obtain the phase correlation array */
fftw_execute(ifft_res);
//normalize and copy to result image
for( i = 0 ; i < fft_size ; i++ ) {
poc_data[i] = res[i][0] / ( double )fft_size;
}
/* deallocate FFTW arrays and plans */
fftw_destroy_plan( fft_img1 );
fftw_destroy_plan( fft_img2 );
fftw_destroy_plan( ifft_res );
fftw_free( img1 );
fftw_free( img2 );
fftw_free( res );
}
Peak FFTW_test(IplImage* src,IplImage* temp)
{
clock_t start=clock();
int t_w=temp->width;
int t_h=temp->height;
/* create a new image, to store phase correlation result */
IplImage* poc = cvCreateImage( cvSize(temp->width,temp->height ), IPL_DEPTH_64F, 1 );
/* get phase correlation of input images */
phase_correlation2D( src, temp, poc );
/* find the maximum value and its location */
CvPoint minloc, maxloc;
double minval, maxval;
cvMinMaxLoc( poc, &minval, &maxval, &minloc, &maxloc, 0 );
/* IplImage* poc_8= cvCreateImage( cvSize(temp->width, temp->height ), 8, 1 );
cvConvertScale(poc,poc_8,(double)255/(maxval-minval),(double)(-minval)*255/(maxval-minval));
cvSaveImage("poc.png",poc_8); */
cvReleaseImage( &poc );
clock_t end=clock();
int time= end-start;
//fprintf( stdout, "Time= %d using clock() \n" ,time/*dt*/ );
//fprintf( stdout, "Maxval at (%d, %d) = %2.4f\n", maxloc.x, maxloc.y, maxval );
CvPoint pt;
pt.x= maxloc.x;
pt.y= maxloc.y;
//4 variants?
//if(maxloc.x>=0&&maxloc.x<=t_w/2&&maxloc.y>=0&&maxloc.y<=t_h/2)
//{
// pt.x= src->width-maxloc.x;
// pt.y= -maxloc.y;
//}
//if(maxloc.x>=t_w/2&&maxloc.x<=t_w&&maxloc.y>=0&&maxloc.y<=t_h/2)
//{
// pt.x= src->width-maxloc.x;
// pt.y= src->height-maxloc.y;
//}
//if(maxloc.x>=0&&maxloc.x<=t_w/2&&maxloc.y>=t_h/2&&maxloc.y<=t_h)
//{
// /*pt.x= -maxloc.x;
// pt.y= -maxloc.y;*/
// pt.x= src->width-maxloc.x;
// pt.y= src->height-maxloc.y;
//}
//if(maxloc.x>=t_w/2&&maxloc.x<=t_w&&maxloc.y>=t_h/2&&maxloc.y<=t_h)
//{
// pt.x= -maxloc.x;
// pt.y= src->height-maxloc.y;
//}
Peak pk;
pk.maxval= maxval;
pk.pt=pt;
return pk;
}