I'm trying to stitch 2 images using cross correlation (phase correlation).Images are the same size.Only shift is present. I tried opencv with cvDFT and opencv+FFTW it seems to work, but for some reason coorelation peak coordinate is not the shift coordinate, maybe it depends on quadrant where correlation point is. So the question is how to obtain shift point from correlation point.
Note: I don't need "interest points" SIFT-like approaches.
My suggestion commented in code. Here is the code that I use:
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;
}
