相位相關算法

相位相關算法：

 1.相位相關簡介：相位相關算法的理論基礎是傅裡葉變換，目前在傅裡葉變換領域有了快速算法fft，比較成熟的庫有fftw開源庫，是以相位相關法有極大的速度優勢，相位相關在圖像融合、模式識别特征比對等有着廣泛應用。

    下面我就圖像融合裡的應用做個簡要介紹：

 針對有平移失配、旋轉的圖像融合分别作介紹。

   1）圖像間有平移變換。

        圖像f2(x,y)是圖像f1(x,y)經平移（x0,y0）後得到的圖像，即

     f2(x,y)=f1(x-x0,y-y0),由傅裡葉時移性質對應傅裡葉變換F1和F2的關系如下：                  F2（u,v）=exp(-j*2*pi(u*x0+v*y0))*F1(u,v)         計算頻域交叉功率譜可得：exp(j*2*pi(u*x0+v*y0))=F1(u,v)*F3 / |F1(u,v)*F3|  F3是F2的共轭。      最後在對交叉功率譜ifft變換可得到一個沖擊函數，此函數在其他位置幾乎為零，隻有在（x0,y0）處有最大值，     是以，可計算出平移參數。   2）針對圖像間有平移旋轉變換關系：      若圖像f2(x,y)是圖像f1(x,y)經平移（x0,y0）、旋轉a角度後得到的圖像，用下面公式表示為：                  f2(x,y)=f1(x*cos(a)+y*sin(a)-x0,-x*sin(a)+y*cos(a)-y0))            由傅裡葉旋轉平移特性，fft變換後兩圖像間的關系如下：                  F2(u,v)=exp(-j2pi(u*x0+v*y0))*F1(u*cos(a)+v*sin(a),-u*sin(a)+v*cos(a))       用M1、M2分别表示F1、F2的能量，則：                     M2（u,v）=M1(u*cos(a)+v*sin(a),-u*sin(a)+v*cos(a));                                 由上式看出F1、F2能量是相同的。把直角坐标轉到極坐标可表示如下：                                   M1（r,a）=M2(r,a-a0)       再由1）所述方法，在極坐标系下用相位相關可求出旋轉角度a0,最後對圖像以角度a0做旋轉，旋轉得到圖像與原圖再次      相位相關就可求出圖像間的平移參數。      下面是對有平移圖像的相位相關的代碼（fft用FFTW庫）：                                                  void phase_correlation( IplImage *ref, IplImage *tpl, IplImage *poc )

{

 int  i, j, k;

 double tmp;  

 int width    = ref->width;

 int height   = ref->height;

 int step     = ref->widthStep;

 int fft_size = width * height;  

 uchar  *ref_data = ( uchar* ) ref->imageData;

 uchar  *tpl_data = ( uchar* ) tpl->imageData;

 double  *poc_data = ( double* )poc->imageData;  

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

 fftw_plan fft_img1 = fftw_plan_dft_1d( width * height, img1, img1, FFTW_FORWARD,  FFTW_ESTIMATE );

 fftw_plan fft_img2 = fftw_plan_dft_1d( width * height, img2, img2, FFTW_FORWARD,  FFTW_ESTIMATE );

 fftw_plan ifft_res = fftw_plan_dft_1d( width * height, res,  res,  FFTW_BACKWARD, FFTW_ESTIMATE );  

 for( i = 0, k = 0 ; i < height ; i++ ) {

  for( j = 0 ; j < width ; j++, k++ ) {

   img1[k][0] = ( double )ref_data[i * step + j];

   img1[k][1] = 0.0;    img2[k][0] = ( double )tpl_data[i * step + j];

   img2[k][1] = 0.0;

  }

 }  

 fftw_execute( fft_img1 );  

 fftw_execute( fft_img2 );  

 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;

 }  

 fftw_execute(ifft_res);  

 for( i = 0 ; i < fft_size ; i++ ) {

  poc_data[i] = res[i][0] / ( double )fft_size;

 }  

 fftw_destroy_plan( fft_img1 );

 fftw_destroy_plan( fft_img2 );

 fftw_destroy_plan( ifft_res );

 fftw_free( img1 );

 fftw_free( img2 );

 fftw_free( res ); 

}