phase correlation using FFTW

Question

Now for phase correlation(two images) I use 1d transform. How to use 2d transform(it will be faster?), how to use wisdom and multythread support? if you give code example is will be better.

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;

    //fftw_init_threads(); // for MT
    //fftw_plan_with_nthreads(2);

    /* 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( width ,height, img1, img1, FFTW_FORWARD,  FFTW_ESTIMATE );
    fftw_plan fft_img2 = fftw_plan_dft_2d( width ,height, img2, img2, FFTW_FORWARD,  FFTW_ESTIMATE );
    fftw_plan ifft_res = fftw_plan_dft_2d( width ,height, res,  res,  FFTW_BACKWARD, FFTW_ESTIMATE );

    /*int f= fftw_init_threads();
    fftw_plan_with_nthreads(10);*/

    /* 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;
        }
    }

    /* 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_cleanup_threads();
    fftw_destroy_plan( fft_img1 );
    fftw_destroy_plan( fft_img2 );
    fftw_destroy_plan( ifft_res );
    fftw_free( img1 );
    fftw_free( img2 );
    fftw_free( res );
}

Answer 1

Finally here is the code:

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;
}  

Answer 2

You will need to use 2D FFTs if you want to calculate phase correlation for two images. You don't need to worry about using FFTW's wisdom for now - just use a basic 2D plan for the FFTs until you get this working. Ditto for multi-threading - get it working single-threaded first.