相位相关算法的详细介绍（一）

相位相关算法：

 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;

 /* get image properties */
 int width    = ref->width;
 int height   = ref->height;
 int step     = ref->widthStep;
 int fft_size = width * height;

 /* setup pointers to images */
 uchar  *ref_data = ( uchar* ) ref->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_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 );

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