相位相关算法的详细介绍(一)

时间:2022-03-10 05:57:46

相位相关算法:

 1.相位相关简介:相位相关算法的理论基础是傅里叶变换,目前在傅里叶变换领域有了快速算法fft,比较成熟的库有fftw开源库,因此相位相关法有极大的速度优势,相位相关在图像融合、模式识别特征匹配等有着广泛应用。

    下面我就图像融合里的应用做个简要介绍:

 针对有平移失配、旋转的图像融合分别作介绍。

   1)图像间有平移变换。

        图像f2(x,y)是图像f1(x,y)经平移(x0,y0)后得到的图像,即

     f2(x,y)=f1(x-x0,y-y0),由傅里叶时移性质对应傅里叶变换F1F2的关系如下:
                 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 ); 
}