一种傅里叶—梅林变换空间图像快速配准算法

焦继超，赵保军，周刚

(1.北京理工大学 信息与电子学院，北京100081;2.中国科学院 长春光学精密机械与物理研究所，吉林 长春130033)

0 引言

为了检测跟踪空间中的暗弱碎片，需要口径大而焦距短的望远镜。但是，随着望远镜口径的扩大，望远镜的焦距也要变长，于是望远镜的CCD 视场就会变小［1］，这就降低了地基天文望远镜探测空间碎片的能力。同时，研制大口径的天文望远镜也存在技术、研制周期等方面的困难。此时，采用并行排列的望远镜阵列能够解决上述问题［2］。由于采用阵列望远镜在同一时刻拍摄的图像并不完全重合，需要对不同CCD 相机拍摄的图像进行配准，才能实现图像的精确叠加，提高图像信噪比，以检测出暗弱的空间碎片。

基于区域的方法是从参考图像中选取一定大小的区域在待配准图像中搜索与其相关性最大的图像块，把2 个区域图像的中心作为一对特征点，然后利用这些特征点求解变换方程［3］。利用灰度值作为相似度测量的算法结构简单［4］，易于实现，但是运算量大，容易受光照、噪声等因素的干扰。采用互信息配准的方法通过计算图像概率密度的相关性实现图像的精确配准，多用于医学图像处理中。但是该算法需要估计图像的概率密度，因此对于概率分布比较复杂的图像具有较大的配准误差。同时，空间图像中的星体基本都呈斑状，结构特征较少，采用基于结构特征的配准方法存在较大困难。

针对一般算法难以实现大分辨率、低信噪比空间图像实时和精确配准的问题，本文采用了一种基于改进的傅里叶-梅林变换的配准算法。实验结果表明本算法能够实现空间图像快速配准的前提下，获得亚像素的配准精度。

1 基于FMT 的图像配准基本原理

基于FMT 的配准算法是从基于FFT 的配准算法发展起来的，它利用图像频谱的相位相关性实现存在平移、旋转和缩放的2 幅图像间的配准［5］。

对于存在平移、旋转和缩放的图像f1(x，y)和f2(x，y)，两者在空域中的位置关系如下所示

其中，f1为参考图像，f2为待配准图像，为了便于叙述，图像在水平和垂直方向采用相同的缩放系数a，θ0为图像的旋转系数，Δx 和Δy 分别为水平方向和垂直方向上的平移量。

对式(1)进行傅里叶变换［6］，得到如下等式

其中:F1(ξ，η)为f1的频谱，F2(ξ，η)为f2的频谱。由上式可知，图像在傅里叶变换前后的旋转角度不变，缩放系数变为原值的倒数。不考虑图像间的平移运动时，F1和F2的幅度值间的关系如下

由上式可知，方程具有三角函数的形式，因此可以将式(3)中的变量进行如下形式的代换

其中，rp(θ，ρ)和sp(θ，ρ)分别为F1和F2在极坐标中的频谱。需要注意的是，频谱的幅度值对应的旋转角度是周期性的，其周期为π，则rp和sp有如下性质［7］

其中，n=…-2，-1，0，1，2，….则式(3)中的坐标变量进行极坐标转换，其形式如下所示

其中，θ 为极坐标系下的角度参数，θ0为图像间的旋转角度。把式(8)和式(9)带入(3)式可得

则上式可以简化为

由上式可知，在极坐标系下，图像旋转角度θ 沿极点进行旋转，缩放系数在极轴方向上成比例变化。因此，对变量进行对数变换

其中:λ 为极坐标系下的极径ρ 的对数，b 为直角坐标系下缩放系数a 的对数。同时我们定义:rpl(θ，λ)=rp(θ，ρ)，spl(θ，λ)=sp(θ，ρ)，则式(10)可以变为如下形式

由式(13)可知，在对数极坐标系下，不但图像旋转角度θ 沿极点进行旋转，而且缩放系数a 也沿极轴进行平移。这2 个性质被称为角度不变性和距离不变性。基于相位相关的方法获得旋转系数和缩放系数，则式(13)两端进行傅里叶变换，

把上式带入交叉能量谱公式中［8］，

当等式左边有最大值时，求得的θ0和b 即为在对数极坐标系中对应的旋转系数和缩放系数。因此可以根据式(16)求得在直角坐标系下的缩放系数

2 新的基于FMT 的空间图像配准算法的原理及实现

2.1 一种新的旋转系数判断方法

空间图像的分辨率一般在1 024 ×1 024 像素或2 048 ×2 048 像素甚至更高，即使采用快速傅里叶变换也要进行大量的计算，影响了空间图像配准的实时性。同时，从第二部分对基于FMT 配准算法原理的论述来看，通过提高图像傅里叶变换的速度来保持配准算法的实时性比较困难，因为目前还没有比FFT 变换更为简单快速的傅里叶变换算法。一方面，拍摄得到的空间图像间存在的平移运动和缩放对空间图像的相关性影响相对较小;另一方面，在频域中，由于旋转系数存在周期性，为了获得正确的系数，需要利用获得的缩放系数和多个旋转系数对待配准图像进行转换并计算其与参考图像的相关系数，这一过程需要大量的计算。

为了提高空间图像配准的速度，保证其实时性，本文提出了一种新的旋转系数判断方法，即通过计算两幅原始配准图像间的相关值T，并与阈值βth进行比较，以获得正确的旋转系数，其具体的判断步骤如下:

1)根据相位相关性，得到在对数极坐标系下配准图像间的旋转系数θ1和θ2且θ1＞θ2;

2)计算参考空间图像和待配准空间图像间的相关值T，公式如下所示［9］，

其中，M 为图像列向量的个数，N 为图像行向量的个数，fs(x，y)为待配准图像，fr(x，y)为参考图像;

3)比较T 与阈值βth间的大小，如果T ＞βth，则说明两幅图像间的相关性较大，选取θ2为旋转系数;如果T＜βth，则说明两幅图像间的差别比较大，选取θ1为旋转系数，其中βth为常数，由大量的实验得到，本文试验中βth取值范围是3 500～3 900.

从以上的步骤来看，本文提出的方法只进行2次相关运算，而原旋转角的判断方法至少要进行4次乘法运算(包括对图像放大σ-1倍时需要插值的2 次运算及计算交叉能量谱的值)和计算交叉能量谱时对图像进行的2 次FFT 变换。以1 024 ×1 024像素的空间图像为例，2 种方法中主要算法的运算量如表1所示。

表1 2 种旋转系数判断方法运算量的比较Tab.1 The comparison of computation between two kinds of rotation coefficient determination methods

表中，表中的除法是计算交叉能量谱时用到的。由表1可以看出，本文提出的判断方法之所以能够减少运算量，增加配准算法的运算速度，就是因为大大减少了计算结构复杂的除法和复数乘法的次数。

2.2 空间图像配准算法中的插值

本文提出的空间图像配准算法需要在不同阶段进行多次插值，采用相同的插值方法会影响整个配准算法的运行速度，因此本文在配准的不同阶段采用不同的插值方法，在不降低配准速度的前提下实现较高的配准精度。

在基于FMT 空间图像配准算法中，为了计算旋转系数和缩放系数，需要进行极坐标和对数坐标的转换，极坐标系和对数极坐标系的坐标点都是非均匀分布的，而直角坐标系坐标点是均匀分布的，直角坐标系中的所有点在转换后不一定都对应其它两种坐标系中的整数点，因此需要进行插值。因为在转换的过程中需要进行2 次插值，为了保证空间图像配准的实时性，需要采用易于实现的插值方法，本文采用的是双线性插值法。

一般情况下，计算得到的转换参数是非整数，当利用转换方程对待配准图像进行转换时，输出像素通常被映射到新坐标系下的非整数位置，因此为了决定与该位置相对应的灰度级，也要进行插值。这个阶段的插值对最后的配准精度有着较大影响，所以采用3 次样条插值。常用的插值方法主要有最邻近插值、双线性插值、样条插值。最邻近插值和双线性插值简单易于实现，但是人工处理的痕迹比较明显，3 次样条插值能够较好的消除锯齿现象，插值精度比较高［10］，3 次样条插值的具体步骤如下，

1)输入m 个插值结点，α =x1＜x2＜…＜xm=β，对应函数值为y1，y2，…，ym，边界条件y'1，y'2，…，y'm，待求插值点x0;

2)计算gi=xi+1-xi，其中i=1，2，…，m-1;

其中，i=2，3，…，m-1;

5)在满足边界条件的前提下，计算如下方程

其中，φn和ωn为插值基函数，

6)输出各区间的3 次样条插值函数si(x).

2.3 空间图像配准的流程图

综合前面的论述，图1给出了本文提出的算法对参考空间图像s(x，y)和待配准空间图像r(x，y)进行配准的流程图，其中θ1＞θ2，如下所示，

图1 空间图像配准流程图Fig.1 The space image registration flowchart

3 实验结果与分析结论

在所研制的DSP+FPGA 硬件平台上通过了实验验证，结构框图如图2所示，其中，实地拍摄的空间图像在 FPGA 中进行预处理，在 DSP(TMS320C6455)中进行基于FMT 的图像配准，DPRAM 存储原始空间图像，SSRAM 缓存图像数据，处理结果发送给主控计算机。

图2 DSP+FPGA 结构框图Fig.2 Block Diagram of DSP+FPGA

为了验证本算法的性能，本文采用实地拍摄且分辨率为1 024 ×1 024 像素的空间图像为参考图像r(x，y)，对r(x，y)进行向下和向右平移(Δx =2，Δy=1)，并作3°的逆时针旋转得到待配准图像s(x，y)，如图3所示，2 幅图像的噪声比较严重，其中标记出的是待配准图像的独有星体。

为了证明本算法在空间图像配准中的性能，不但采用均方根误差(RMSE)准则对实验结果进行评价，而且还同基于区域配准和基于结构特征的配准算法进行比较。

图3 配准空间图像Fig.3 The registration space images

3.1 本文配准算法性能分析

由图3可以发现，空间图像的信噪比很低，有着很多明显的噪点，但是通过高通滤波可增强星体边缘，抑制部分噪声，图4给出了图3经过傅里叶变换和滤波后的三维频谱图。

图4 空间图像直角坐标系频谱图Fig.4 The space image spectrum in cartesian coordinate

从图4可以看出，空间图像频谱中某一个频率对应一个全局的峰值，而噪声对应的幅度值要小于这个值，因此可利用2 个频谱幅度值峰值的相关性实现图像的配准。

为了获得配准转换方程的旋转系数和缩放系数，需要将直角坐标系下的频谱转换到对数极坐标系下，图5给出了对数极坐标系下的空间图像频谱图。

由图5可知，对数极坐标系下的空间图像频谱和直角坐标系下的频谱类似，某一坐标值对应一个全局最大幅度值。

利用相位相关性，分别在对数极坐标系和直角坐标系中求出旋转角度θ =3.02°、缩放系数a =1.01 和平移系数(Δx =2.00，Δy =1.00)，图6给出了配准后的图像s'(x，y)以及和r(x，y)的叠加图像c(x，y).

从叠加图像c(x，y)中可以看出，不论是面积很小的星体还是拖尾都能够完全重合，而没有明显的叠加痕迹，其中圆形标记的和图3(b)中标记出的是同一星体。

同时本文还采用均方根误差(RMSE)准则来评价算法性能，其定义如下［11］

图5 空间图像对数极坐标系频谱图Fig.5 The space image spectrum in log-polar coordinate

其中，N 为空间图像提取的特征星体个数，(x，y)为参考图像中特征星体的质心位置，(x'，y')为配准后图像中特征星体的质心位置。

利用式(19)可以求得空间图像的RMSE，配准前，RMSE=1.202，配准后，RMSE=0.518，即配准前后RMSE 减少值ΔRMSE =0.684.从主观观测和客观指标方面都说明了本文配准算法对空间图像进行了精确配准。

3.2 本文算法与其它配准算法的性能比

为了进一步突出本算法的性能，本次试验又利用基于区域配准算法和基于结构特征提取配准算法进行空间图像的配准，表2在运算时间和RMSE 两个指标上对3 种配准算法进行了比较。

由表2可知，在运算速度上，本文提出的算法采用了新的旋转角度判断方法和多种插值方式，因此分别比基于区域和基于结构特征的配准算法分别快4.155 s 和3.476 s;在RMSE 上，基于FMT 的配准算法比其他两种配准算法分别减少0.436 s 和0.188 s.

图6 配准后和叠加空间图像Fig.6 The space image after registration and superposable image

表2 3 种配准算法的性能比较Tab.2 The comparison of the three registration methods

4 结论

本文根据空间图像配准实时性和精确性的要求，提出了一种新的基于FMT 的配准算法。首先，提出一种新的基于原始空间图像相关性的旋转系数判断方法。然后，为了在不损失配准精度前提下降低算法复杂度，提出在坐标系转换和待配准图像变换阶段采用不同插值方法，并给出了三次样条插值的步骤。最后，给出了本文算法实现空间图像配准的流程图。通过硬件平台的实验验证，结果证明:基于FMT 的配准方法能够基本实现空间图像配准的实时性要求;配准精度达到0.5 像素;算法的RMSE减少0.684.基本满足了空间图像配准对实时性和精确性的要求。

References)

［1］王鸣浩，陈涛，王建立，等.捆绑式望远镜图像信噪比测量及分析［J］.光学精密工程，2009，17(1):92 -97.WANG Ming-hao，CHEN Tao，WANG Jian-li，et al.Measurement and analysis of image SNR in binding style telescope［J］.Opt.Precision Eng.，2009，17(1):92 -97.(in Chinese)

［2］Nicholas K.Pan-Starrs:a wide-field optical survey telescope array［J］.SPIE，2004，5489:230 -248.

［3］Héctor Fernando Gómez-García，José L.Marroquín a，Johan Van Horebeek.Image registration based on kernel-predictability［J］.Computer Vision and Image Understanding，2008，112:160 -172.

［4］Alliney S，Cortelazzo G，Mian G.A.On the registrations of an object translating on a static background ［J］.Pattern Recognition，1996，29(1):131 -141.

［5］刘卫光，崔江涛，周利华.插值和相位相关的图像亚像素配准方法［J］.计算机辅助设计与图形学学报，2005，17(6):1273-1277.LIU Wei-guang，CUI Jiang-tao，ZHOU Li-hua.Subpixel registration based on interpolation and extension phase correlation［J］.Journal of Computer-Aided Design ＆ Computer Graphics，2005，17(6):1273 -1277.(in Chinese)

［6］高莹莹，杨建峰，马晓龙，等.基于Fourier-Mellin 算法的干涉图像配准［J］.光学精密工程，2007，15(9):1415 -1420.GAO Ying-ying，YANG Jian-feng，MA Xiao-long，et al.Interference image registration based on Fourier-Mellin algorithm［J］.Opt.Precision Eng.，2007，15(9):1415 -1420.(in Chiese)

［7］Alexander W，Paul F.Fast phase-based registration of multimodal image data［J］.Signal Processing，2009，89:724 -737.

［8］Barbara Zitová，Jan Flusser.Image registration methods:a survey［J］.Image and Vision Computing，2003，21:977 -1000.

［9］Kamp S，Heyden D，Ohm J R.Inter-temporal vector prediction for motion estimation in scalable video coding［C］∥Intelligent Signal Processing and Communication Systems，ISPACS 2007.Xiamen，2007:586 -589.

［10］Wang Q，Ward R K.A new orientation-adaptive interpolation method［J］.IEEE Transactions on Image Processing，2007，16:889 -990.

［11］Chahira Serief，Mourad Barkat，Youcef Bentoutou.Elastic registration of remote-sensing images based on the nonsubsampled contourlet transform［J］.Transactions on Pattern Analysis and Machine Intelligence，IEEE，2007，14(3):1021 -102.