相位相关使用FFTW

现在对于相位相关（两个图像）我使用1D变换。如何使用2D变换（它会更快？），如何使用智慧和多线程支持？如果你给代码示例是会更好。相位相关使用FFTW

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

来源

2012-02-03 mrgloom

最后这里是代码：

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

来源

2012-08-20 09:24:57 mrgloom

请注意，您可以只用3次而不是4次进行复数乘法运算。尽管如此，我认为它不会造成很大的速度差异。 http://mathworld.wolfram.com/ComplexMultiplication.html – CookieOfFortune 2012-10-08 16:17:17

如果要计算两个图像的相位相关性，则需要使用2D FFT。您现在不需要担心使用FFTW的智慧 - 只需使用FFT的基本2D计划，直到获得此效果。同样适用于多线程 - 首先让它工作在单线程上。

来源

2012-02-03 09:33:24

我实现它像（http://codepaste.ru/9225/），但似乎它的工作不一样的1D变换。 – mrgloom 2012-02-03 10:29:32

我看不出有什么明显的错误 - 尝试用IMG1测试它== IMG2 - 你应该在（0,0）的单峰。 – 2012-02-03 10:35:14

它给我（0,0），如果我把它叫做phase_correlation2D（src，src，res）;但它在真实测试图像上给出的结果与phase_correlation1D不一样。另外通过2D变换的方式更快。 – mrgloom 2012-02-03 10:57:42

相位相关使用FFTW

回答

相关问题