相位相关算法

相位相关算法：

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

}