一种高精度的亚像素级图像配准算法

周玫君，赵辽英，厉小润

(1.杭州电子科技大学计算机学院，浙江 杭州 310018；2.浙江大学信息与电子工程系，浙江 杭州 310027)

0 引 言

图像配准是影响图像拼接、图像变化检测和图像融合等图像处理的关键技术[1]，广泛应用于医学图像处理、高分辨率图像生成以及图像重建等实际应用中，要求配准精度达到亚像素级[2]。现有的亚像素级图像配准算法主要有基于灰度插值法[3]、解最优化问题法[4]和变换域求解法[5-7]。基于灰度插值法首先通过再采样使2幅图像放大并且具有相同的分辨率，然后进行配准，但在再采样时必须进行插值，而插值算法的质量决定了插值方法的精确度。解最优化问题法可归结为求配准参数最小化的代价函数问题，适用于多光谱图像及亮度变化的配准，但算法复杂度较高。变换域求解法对噪声、光照等鲁棒性强，对图像灰度信息依赖较低，因此受到广泛关注。早期的变换域求解法主要采用扩展相位相关法，如文献[5]利用曲线拟合的思想将相位相关算法推广到亚像素精度，但达到的精度有限并且只能解决平移问题。为求解旋转和缩放参数，常用的方法是相位相关法结合扩展傅里叶梅林变换(Fourier-Mellin，FM)法，如文献[6]提出FM结合sinc函数来近似表示图像间的相位相关函数，在一定程度上提高了运动参数的估计精度，但该方法受限于sinc函数和相位相关函数之间的误差，并且没有考虑缩放问题。为进一步提高算法的精度，袁恒等[7]将矩阵乘法离散傅里叶变换与FM结合提出上采样快速FM变换法，但该方法不能实现任意精度的亚像素参数估计。针对高精度旋转缩放参数估计问题，本文综合考虑配准精度和运算速度，提出一种改进的曲面拟合修正上采样傅里叶梅林算法，在不降低执行效率的情况下，可获得任意高精度的亚像素级运动参数。

1 基础理论

基于傅里叶变换求图像间的旋转、缩放和平移参数主要使用2个技术：相位相关匹配[8]和傅里叶梅林变换[9]。

1.1 相位相关匹配

设大小均为M×N图像f1(x,y)和f2(x,y)存在如下的平移关系：

f2(x,y)=f1(x-x0,y-y0)

(1)

则f1(x,y)和f2(x,y)之间的归一化互功率谱为：

(2)

对式(2)进行傅里叶反变换可得f1(x,y)和f2(x,y)之间的脉冲函数δ：

q(x,y)=δ(x-x0,y-y0)

(3)

搜索其峰值位置，即可确定图像f2相对图像f1的像素级平移参数(x0,y0)。

1.2 Fourier-Mellin变换

假设2幅图像f1(x,y)，f2(x,y)存在平移、旋转和缩放变换：

(4)

式中，φ，s，Δx，Δy分别为旋转、缩放和平移参数，则f1(x,y)，f2(x,y)的傅里叶变换的幅度谱M1，M2满足如下关系：

(5)

将幅度谱转换到对数-极坐标中，可以得到：

M2(ε,θ)=M1(ε-d,θ-φ)

(6)

2 曲面拟合修正上采样傅里叶梅林法

在实际的应用场景中，图像间的运动都是连续的，如果采用第1节中的相位相关法和傅里叶梅林法估计图像运动，只能得到对数-极坐标下的整像素级运动参数，不能解决图像的实际问题。为了得到任意精度的亚像素级运动参数，本文提出改进的曲面拟合修正上采样傅里叶梅林算法，首先，对预处理后的2幅图像进行傅里叶变换，并对其归一化互功率谱Q(u,v)执行2倍零填充，反变换检测其峰值坐标得到旋转缩放矢量的1/2像素级估计值(dε2,dθ2)；然后，应用矩阵乘法离散傅里叶变换快速计算n倍上采样相位相关峰作为旋转缩放矢量的亚像素初始估计值；最后，采用曲面拟合法对亚像素初始估计值进行修正，进而得到任意精度的运动估计结果。

2.1 亚像素级旋转缩放初始估计值的快速计算

(7)

2.2 曲面拟合计算旋转缩放参数调整值

在n倍上采样相位相关曲面上，采用二元二次多项式曲面拟合函数在以初始估计值(ε0,θ0)为中心的3×3领域进行拟合：

q(ε,θ)=a0+a1ε+a2θ+a3ε2+a4εθ+a5θ2

(8)

式中，q(ε,θ)表示坐标为(ε,θ)的相关函数值。将点(ε0,θ0)与其周围8点的相关函数值代入式(8)，并利用最小二乘法得到该多项式的系数a0,a1,a2,a3,a4,a5，即可求出拟合曲面的极值点位置：

(9)

用拟合曲面的极值点坐标(Δε,Δθ)除以n的结果(dεn,dθn)对1/2像素级的估计值(dε2,dθ2)进行更新，得到对数-极坐标下的任意精度的亚像素级旋转、缩放估计值，ε=dε2+dεn,θ=dθ2+dθn。再对ε,θ进行坐标变换公式计算得到旋转参数φ和缩放参数s。

图1 算法总体流程

2.3 算法总体流程

按照上述分析，算法总体流程如图1所示，主要步骤如下：

(1)对图像进行预处理。为解决频谱混叠和傅里叶变换的边缘效应对配准精度的影响，首先对图像进行梯度[10]和加窗处理[11]。

(2)计算初始旋转缩放像素级参数。对预处理后的2幅图像进行傅里叶变换并对其归一化互功率谱执行2倍零填充，反变换检测其峰值坐标即可得到旋转缩放矢量的1/2像素级估计值(dε2,dθ2)。

(3)求n倍上采样相位相关峰。利用矩阵乘法离散傅里叶变换[9]快速计算n倍上采样相位相关曲面，求得的峰值坐标(ε0,θ0)即为坐标反变换前的旋转缩放亚像素级初始估计参数。

(4)曲面拟合计算旋转缩放参数调整值。在上采样相位相关曲面上，对以峰值坐标为中心的3×3区域进行曲面拟合，修正峰值以得到任意精度的亚像素级调整值(dεn,dθn)。

(5)使用(dεn,dθn)对初始估计值(dε2,dθ2)进行修正，从而实现任意精度级别的旋转缩放运动估计。

最后根据对数-极坐标和笛卡尔坐标间的映射，计算笛卡尔坐标下的旋转缩放参数，对待配准图进行反变换，再利用基于拟合的扩展相位相关法求亚像素平移参数。

3 实验结果与分析

分别采用仿真图像配准实验和真实图像配准实验来验证本文算法的性能。选用FM结合sinc曲线拟合的方法[6](简称FM_sinc)、快速上采样FM变换法[7](简称Fast_FM)和本文提出的曲面拟合修正上采样FM算法进行实验比较。由于本文算法的重点在于旋转和缩放参数的精确估计，因此，实验中所有方法采用相同的算法估计平移参数，预处理步骤一致。实验环境配置为：Inter CPU(i7-7700,3.60 GHz)、8.0 GB内存、64位的Windows10操作系统。

3.1 仿真序列图像实验

图2(a)为某地区上空航拍遥感影像中截取的大小为330×330的子图。以图2(a)为参考图，取随机值(3.6，2.4)平移变换后，依次对其进行步长为5°的旋转操作，缩放系数每隔2幅图增加0.1，共得到6幅序列图像，图2(b)和(c)分别为第3和第6幅序列图。

图2 仿真实验的部分图像

用配准后的均方根误差(Root Mean Square Error，RMSE)评价配准精度，其表达式为：

(10)

图3 3种算法的配准均方根误差

分别采用3种算法求得的配准后图像的RRMSE如图3所示。经计算，FM_sinc，Fast_FM和本文算法的平均RRMSE分别为0.862 0，0.427 3和0.154 1，本文算法的配准精度分别提高了0.707 9像素和0.273 2像素。从图3可以看出：FM_sinc算法的RRMSE最大，并随着图像旋转缩放的增大而增大。Fast_FM算法对旋转缩放参数的估计精度有明显的提升，但其精度受限于n的取值，只能得到1/n精度估计结果。而本文算法在1/n精度估计的基础上，利用拟合修正估计值，实现了任意精度估计结果。

实验中，通过计算可得，FM_sinc，Fast_FM和本文算法的平均运行时间分别为0.763 7 s，0.790 1 s和0.805 3 s，算法运行时间分别增加约0.041 6 s和0.015 2 s。说明本文算法能在增加较少运行时间的情况下获得更高的配准精度。

仿真图像在平移固定，旋转缩放(15°，1.2)下，采用3种算法配准后的图像与参考图像之间的联合概率直方图如图4所示。

图4 旋转缩放(15°，1.2)时，3种算法配准结果的联合直方图

从图4可以看出：图4(a)中有大量的点分散在对角线周围，相比图4(a)，图4(b)中有部分点更靠近对角线，图4(c)中点的分布最靠近直角线。联合直方图中的点分布越靠近对角线，误差越小，图4结果表明本文算法的配准误差最小。

3.2 真实图像实验

图5 真实图像

选取某地区2001年拍摄的SPOT4图像和2004年拍摄的Landsat TM图像作为真实图像进行实验。分别从中心截取大小为330×330的子图作为参考图和待配准图，如图5所示。

由于真实遥感图像间的旋转、缩放和平移量未知，无法对结果客观评价，所以根据配准后的拼接效果进行主观评价。图6给出3种算法配准后的拼接结果图，其中每幅图右边的小图依次为左图中从左到右的圆圈中拼接区域的局部放大图。从图6可以看出，3种算法都得到了很好的配准结果，但仔细分析各图的第一个局部放大区域可以看出，FM_sinc求得的配准后图像在拼接处有完全错开的断层，Fast_FM算法求得的图像拼接处断层错开近1/4，本文算法比Fast_FM的拼接效果稍微好一些，断层处错开了近1/2。由此可见，本文算法的配准效果最佳。

图6 3种算法配准拼接结果

4 结束语

针对亚像素级的旋转、缩放运动的高精度估计问题，本文提出一种基于曲面拟合修正矩阵乘法的傅里叶梅林变换算法，计算结果表明本文算法的计算精度误差为0.15像素级别，是一种计算精度高的配准方法，具有广泛的应用前景。如何进一步提高计算精度，增强其鲁棒性，减少计算时间是下一步研究的重点。