数字图像处理，相位相关图像配准算法的C++实现

参考论文中的文字：图像配准是图像处理的基本任务之一,用于将不同时间、不同传感器、不同视角及不同拍摄条件下获取的关于同一目标或场景的两幅或多幅图像进行主要是几何意义上的匹配套和的过程。在对图像配准的研究过程中,大量技术被应用于针对不同数据和问题的图像配准工作,产生了多种不同形式的图像配准技术。

图像配准的基本问题是找出一种图像转换方法,用以纠正图像的形变。造成图像形变的原因多种多样,例如对于遥感图像而言,传感器噪声、由传感器视点变化或平台不稳定造成的透视变化、被拍摄物体的移动、变形或生长等变化、闪电和大气变化,以及阴影和云层遮盖都使图像产生不同形式的形变。正是图像形变原因和形式的不同决定了多种多样的图像配准技术。

迄今已报道了多种图像配准方法,主要有互相关法、傅立叶变换法、点映射法口脚外和弹性模型法。其中傅立叶变换法基于傅立叶变换的相位匹配是利用傅立叶变换的性质而出现的一种图像配准方法。图像经过傅立叶变换,由空域变换到频率缘则两组数据在空何上的相关运算可以变为频谱的复数乘法运算,同时图像在变换域中还能获得在空域中很难获得的特征。

一，基于相位相关的图像配准方法

在时域中信号的平移运动可以通过在频域中相位的变化表现出来（这是傅里叶变换的特性）。同理,图像的旋转、平移和比例变化也能在傅里叶变换的频域中反映出来。而且使用频域方法的好处是计算简单,同时傅立叶变换可以采用方法提高执行的速度。因此,傅氏变换是图像配准中常用的方法之一。下面我们就具体分析当图像发生平移、旋转和缩放时,图像信号在频域中的表现。

1，原理阐述

通过求取互功率谱的傅立叶反变换,得到一个狄拉克函数（脉冲函数）,再寻找函数峰值点对应的坐标,即可得到我们所要求得的配准点。实际上,在计算机处理中,连续域要用离散域代替,这使得狄拉克函数转化为离散时间单位冲击函数序列的形式。在实际运算中,两幅图像互功率谱相位的反变换,总是含有一个相关峰值代表两幅图像的配准点,和一些非相关峰值,相关峰值直接反映两幅图像间的一致程度。更精确的讲,相关峰的能量对应重叠区域的所占百分比,非相关峰对应非重叠区域所占百分比。由此我们可以看出,当两幅图像重叠区域较小时,采用本方法就不能检测出两幅图像的平移量。

2，配准情况

当图像间仅存在平移时，正确的配准图像如图a所示（中心平移化了），最大峰的位置就是两图像的相对平移量，反之若不存在单纯的平移，则会出现如b所示的情况（多脉冲林立）

3，算法流程

二，程序实现

1，核心函数的实现：

[html]  view plain copy print ?

1. void PhaseCorrelation2D(const BYTE *signal,//原信号  
2.     const BYTE *pattern,//带配准信号  
3.     int &height_offset,//高的偏移量  
4.     int &width_offset)//宽的偏移量  
5. {  
6.   
7.     fftw_complex *signal_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);  
8.     fftw_complex *pattern_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);  
9.   
10.     for (int i = 0; i < nRow*nCol; i++)  
11.     {  
12.         signal_img[i][0] = signal[i];  
13.         signal_img[i][1] = 0;  
14.     }  
15.     for (int j = 0; j < nRow*nCol; j++)  
16.     {  
17.         pattern_img[j][0] = pattern[j];  
18.         pattern_img[j][1] = 0;  
19.     }  
20.   
21.     // 对两幅图像傅里叶变换  
22.     fftw_plan signal_forward_plan = fftw_plan_dft_2d(nRow, nCol, signal_img, signal_img,  
23.         FFTW_FORWARD, FFTW_ESTIMATE);  
24.     fftw_plan pattern_forward_plan = fftw_plan_dft_2d(nRow, nCol, pattern_img, pattern_img,  
25.         FFTW_FORWARD, FFTW_ESTIMATE);  
26.     fftw_execute(signal_forward_plan);  
27.     fftw_execute(pattern_forward_plan);  
28.   
29.     // 求互功率谱  
30.     fftw_complex *cross_img = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nRow*nCol);  
31.     double temp;  
32.     for (int i = 0; i < nRow*nCol; i++)  
33.     {  
34.         cross_img[i][0] = (signal_img[i][0] * pattern_img[i][0]) -   
35.             (signal_img[i][1] * (-1.0*pattern_img[i][1]));  
36.         cross_img[i][1] = (signal_img[i][0] * (-1.0*pattern_img[i][1])) +   
37.             (signal_img[i][1] * pattern_img[i][0]);  
38.         temp = sqrt(cross_img[i][0] * cross_img[i][0] + cross_img[i][1] * cross_img[i][1]) + 0.001;  
39.         cross_img[i][0] /= temp;  
40.         cross_img[i][1] /= temp;  
41.     }  
42.   
43.     //对互功率谱求反变换  
44.     fftw_plan cross_backward_plan = fftw_plan_dft_2d(nRow, nCol, cross_img, cross_img,  
45.         FFTW_BACKWARD, FFTW_ESTIMATE);  
46.     fftw_execute(cross_backward_plan);  
47.   
48.     // 释放内存  
49.     fftw_destroy_plan(signal_forward_plan);  
50.     fftw_destroy_plan(pattern_forward_plan);  
51.     fftw_destroy_plan(cross_backward_plan);  
52.     fftw_free(signal_img);  
53.     fftw_free(pattern_img);  
54.   
55.     double *cross_real=new double[nRow*nCol];  
56.     for (int i = 0; i < nRow*nCol; i++)  
57.         cross_real[i] = cross_img[i][0];  
58.   
59.     int max_loc = 0;//准备存放最大值的位置坐标（注意，只有一个值）  
60.     double max_vlaue = 0.0;  
61.     for (int i = 0; i < nRow*nCol; i++)  
62.     {  
63.         if (max_vlaue<cross_real[i])  
64.         {  
65.             max_vlaue = cross_real[i];  
66.             max_loc = i;  
67.         }  
68.     }  
69.   
70.     height_offset = floor(((int)max_loc) / nCol);  
71.     width_offset = (int)max_loc - nCol*height_offset;  
72.   
73.     if (height_offset > 0.5*nRow)  
74.         height_offset = height_offset - nRow;  
75.     if (width_offset > 0.5*nCol)  
76.         width_offset = width_offset - nCol;  
77.   
78.     delete[] cross_real;  
79.     cross_real = NULL;  
80. }  

2，配准模拟结果：

3，实际的图像配准效果

注：下面这张图是单帧图像效果，有大量的椒盐噪声和散斑噪声

注：下面这张图像是八张单帧图像（图像内容可能一样，但是可能存在较小的位移）的叠加效果，散斑噪声和椒盐噪声得到明显去除

注：下面这张图像是八张图像先配准再叠加平均的效果，背景噪声及其散斑噪声得到消除，并且图像对比度得到一定程度提高

注：本博文为EbowTang原创，后续可能继续更新本文。如果转载，请务必复制本条信息！

原文地址：http://blog.csdn.net/ebowtang/article/details/47816245

原作者博客：http://blog.csdn.net/ebowtang