template matching (фазовая корреляция)

[mrgloom_]

Склеиваю два изображения которые имеют некоторое перекрытие (допустим перекрываются на 20%) использую фазовую корреляцию.
точку максимальной корреляции высчитываю по формуле.
$pos(x,y)= Max[Norm[IFFT[FFT[img1]*FFT[img2]]]] $
изображения одного размера, помогите интерпретировать находимый пик корреляции, т.е. из пика корреляции я хочу получить точку смещения.

попытался сравнить код через старый opencv и код через FFTW+opencv
похоже что оба работают правильно, только постоянно какие то непонятные доп. смещения по краям(в коде они закомментированы прибавляются или вычитаются от конечной точки), т.е. я так понимаю в зависимости от того в какой квадрант попадает пик корреляции зависит какая будет формула расчета конечного результата, причем для 2-х методов формулы вроде бы будут разные.
эти коды я использую не для поиска темплейта, а для сшивки пары изображений.изображения одного размера.

вообщем вопрос стоит так как по точке выдаваемой из функции корреляции получить точку смещения?

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 случая локации точки максимальной кореляции
   //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;
}