一種高精度的亞像素級圖像配準(zhǔn)算法

周玫君，趙遼英，厲小潤

(1.杭州電子科技大學(xué)計算機(jī)學(xué)院，浙江 杭州 310018；2.浙江大學(xué)信息與電子工程系，浙江 杭州 310027)

0 引 言

圖像配準(zhǔn)是影響圖像拼接、圖像變化檢測和圖像融合等圖像處理的關(guān)鍵技術(shù)[1]，廣泛應(yīng)用于醫(yī)學(xué)圖像處理、高分辨率圖像生成以及圖像重建等實際應(yīng)用中，要求配準(zhǔn)精度達(dá)到亞像素級[2]?，F(xiàn)有的亞像素級圖像配準(zhǔn)算法主要有基于灰度插值法[3]、解最優(yōu)化問題法[4]和變換域求解法[5-7]?；诨叶炔逯捣ㄊ紫韧ㄟ^再采樣使2幅圖像放大并且具有相同的分辨率，然后進(jìn)行配準(zhǔn)，但在再采樣時必須進(jìn)行插值，而插值算法的質(zhì)量決定了插值方法的精確度。解最優(yōu)化問題法可歸結(jié)為求配準(zhǔn)參數(shù)最小化的代價函數(shù)問題，適用于多光譜圖像及亮度變化的配準(zhǔn)，但算法復(fù)雜度較高。變換域求解法對噪聲、光照等魯棒性強(qiáng)，對圖像灰度信息依賴較低，因此受到廣泛關(guān)注。早期的變換域求解法主要采用擴(kuò)展相位相關(guān)法，如文獻(xiàn)[5]利用曲線擬合的思想將相位相關(guān)算法推廣到亞像素精度，但達(dá)到的精度有限并且只能解決平移問題。為求解旋轉(zhuǎn)和縮放參數(shù)，常用的方法是相位相關(guān)法結(jié)合擴(kuò)展傅里葉梅林變換(Fourier-Mellin，F(xiàn)M)法，如文獻(xiàn)[6]提出FM結(jié)合sinc函數(shù)來近似表示圖像間的相位相關(guān)函數(shù)，在一定程度上提高了運(yùn)動參數(shù)的估計精度，但該方法受限于sinc函數(shù)和相位相關(guān)函數(shù)之間的誤差，并且沒有考慮縮放問題。為進(jìn)一步提高算法的精度，袁恒等[7]將矩陣乘法離散傅里葉變換與FM結(jié)合提出上采樣快速FM變換法，但該方法不能實現(xiàn)任意精度的亞像素參數(shù)估計。針對高精度旋轉(zhuǎn)縮放參數(shù)估計問題，本文綜合考慮配準(zhǔn)精度和運(yùn)算速度，提出一種改進(jìn)的曲面擬合修正上采樣傅里葉梅林算法，在不降低執(zhí)行效率的情況下，可獲得任意高精度的亞像素級運(yùn)動參數(shù)。

1 基礎(chǔ)理論

基于傅里葉變換求圖像間的旋轉(zhuǎn)、縮放和平移參數(shù)主要使用2個技術(shù)：相位相關(guān)匹配[8]和傅里葉梅林變換[9]。

1.1 相位相關(guān)匹配

設(shè)大小均為M×N圖像f1(x,y)和f2(x,y)存在如下的平移關(guān)系：

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

(1)

則f1(x,y)和f2(x,y)之間的歸一化互功率譜為：

(2)

對式(2)進(jìn)行傅里葉反變換可得f1(x,y)和f2(x,y)之間的脈沖函數(shù)δ：

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

(3)

搜索其峰值位置，即可確定圖像f2相對圖像f1的像素級平移參數(shù)(x0,y0)。

1.2 Fourier-Mellin變換

假設(shè)2幅圖像f1(x,y)，f2(x,y)存在平移、旋轉(zhuǎn)和縮放變換：

(4)

式中，φ，s，Δx，Δy分別為旋轉(zhuǎn)、縮放和平移參數(shù)，則f1(x,y)，f2(x,y)的傅里葉變換的幅度譜M1，M2滿足如下關(guān)系：

(5)

將幅度譜轉(zhuǎn)換到對數(shù)-極坐標(biāo)中，可以得到：

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

(6)

2 曲面擬合修正上采樣傅里葉梅林法

在實際的應(yīng)用場景中，圖像間的運(yùn)動都是連續(xù)的，如果采用第1節(jié)中的相位相關(guān)法和傅里葉梅林法估計圖像運(yùn)動，只能得到對數(shù)-極坐標(biāo)下的整像素級運(yùn)動參數(shù)，不能解決圖像的實際問題。為了得到任意精度的亞像素級運(yùn)動參數(shù)，本文提出改進(jìn)的曲面擬合修正上采樣傅里葉梅林算法，首先，對預(yù)處理后的2幅圖像進(jìn)行傅里葉變換，并對其歸一化互功率譜Q(u,v)執(zhí)行2倍零填充，反變換檢測其峰值坐標(biāo)得到旋轉(zhuǎn)縮放矢量的1/2像素級估計值(dε2,dθ2)；然后，應(yīng)用矩陣乘法離散傅里葉變換快速計算n倍上采樣相位相關(guān)峰作為旋轉(zhuǎn)縮放矢量的亞像素初始估計值；最后，采用曲面擬合法對亞像素初始估計值進(jìn)行修正，進(jìn)而得到任意精度的運(yùn)動估計結(jié)果。

2.1 亞像素級旋轉(zhuǎn)縮放初始估計值的快速計算

(7)

2.2 曲面擬合計算旋轉(zhuǎn)縮放參數(shù)調(diào)整值

在n倍上采樣相位相關(guān)曲面上，采用二元二次多項式曲面擬合函數(shù)在以初始估計值(ε0,θ0)為中心的3×3領(lǐng)域進(jìn)行擬合：

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

(8)

式中，q(ε,θ)表示坐標(biāo)為(ε,θ)的相關(guān)函數(shù)值。將點(ε0,θ0)與其周圍8點的相關(guān)函數(shù)值代入式(8)，并利用最小二乘法得到該多項式的系數(shù)a0,a1,a2,a3,a4,a5，即可求出擬合曲面的極值點位置：

(9)

用擬合曲面的極值點坐標(biāo)(Δε,Δθ)除以n的結(jié)果(dεn,dθn)對1/2像素級的估計值(dε2,dθ2)進(jìn)行更新，得到對數(shù)-極坐標(biāo)下的任意精度的亞像素級旋轉(zhuǎn)、縮放估計值，ε=dε2+dεn,θ=dθ2+dθn。再對ε,θ進(jìn)行坐標(biāo)變換公式計算得到旋轉(zhuǎn)參數(shù)φ和縮放參數(shù)s。

圖1 算法總體流程

2.3 算法總體流程

按照上述分析，算法總體流程如圖1所示，主要步驟如下：

(1)對圖像進(jìn)行預(yù)處理。為解決頻譜混疊和傅里葉變換的邊緣效應(yīng)對配準(zhǔn)精度的影響，首先對圖像進(jìn)行梯度[10]和加窗處理[11]。

(2)計算初始旋轉(zhuǎn)縮放像素級參數(shù)。對預(yù)處理后的2幅圖像進(jìn)行傅里葉變換并對其歸一化互功率譜執(zhí)行2倍零填充，反變換檢測其峰值坐標(biāo)即可得到旋轉(zhuǎn)縮放矢量的1/2像素級估計值(dε2,dθ2)。

(3)求n倍上采樣相位相關(guān)峰。利用矩陣乘法離散傅里葉變換[9]快速計算n倍上采樣相位相關(guān)曲面，求得的峰值坐標(biāo)(ε0,θ0)即為坐標(biāo)反變換前的旋轉(zhuǎn)縮放亞像素級初始估計參數(shù)。

(4)曲面擬合計算旋轉(zhuǎn)縮放參數(shù)調(diào)整值。在上采樣相位相關(guān)曲面上，對以峰值坐標(biāo)為中心的3×3區(qū)域進(jìn)行曲面擬合，修正峰值以得到任意精度的亞像素級調(diào)整值(dεn,dθn)。

(5)使用(dεn,dθn)對初始估計值(dε2,dθ2)進(jìn)行修正，從而實現(xiàn)任意精度級別的旋轉(zhuǎn)縮放運(yùn)動估計。

最后根據(jù)對數(shù)-極坐標(biāo)和笛卡爾坐標(biāo)間的映射，計算笛卡爾坐標(biāo)下的旋轉(zhuǎn)縮放參數(shù)，對待配準(zhǔn)圖進(jìn)行反變換，再利用基于擬合的擴(kuò)展相位相關(guān)法求亞像素平移參數(shù)。

3 實驗結(jié)果與分析

分別采用仿真圖像配準(zhǔn)實驗和真實圖像配準(zhǔn)實驗來驗證本文算法的性能。選用FM結(jié)合sinc曲線擬合的方法[6](簡稱FM_sinc)、快速上采樣FM變換法[7](簡稱Fast_FM)和本文提出的曲面擬合修正上采樣FM算法進(jìn)行實驗比較。由于本文算法的重點在于旋轉(zhuǎn)和縮放參數(shù)的精確估計，因此，實驗中所有方法采用相同的算法估計平移參數(shù)，預(yù)處理步驟一致。實驗環(huán)境配置為：Inter CPU(i7-7700,3.60 GHz)、8.0 GB內(nèi)存、64位的Windows10操作系統(tǒng)。

3.1 仿真序列圖像實驗

圖2(a)為某地區(qū)上空航拍遙感影像中截取的大小為330×330的子圖。以圖2(a)為參考圖，取隨機(jī)值(3.6，2.4)平移變換后，依次對其進(jìn)行步長為5°的旋轉(zhuǎn)操作，縮放系數(shù)每隔2幅圖增加0.1，共得到6幅序列圖像，圖2(b)和(c)分別為第3和第6幅序列圖。

圖2 仿真實驗的部分圖像

用配準(zhǔn)后的均方根誤差(Root Mean Square Error，RMSE)評價配準(zhǔn)精度，其表達(dá)式為：

(10)

圖3 3種算法的配準(zhǔn)均方根誤差

分別采用3種算法求得的配準(zhǔn)后圖像的RRMSE如圖3所示。經(jīng)計算，F(xiàn)M_sinc，F(xiàn)ast_FM和本文算法的平均RRMSE分別為0.862 0，0.427 3和0.154 1，本文算法的配準(zhǔn)精度分別提高了0.707 9像素和0.273 2像素。從圖3可以看出：FM_sinc算法的RRMSE最大，并隨著圖像旋轉(zhuǎn)縮放的增大而增大。Fast_FM算法對旋轉(zhuǎn)縮放參數(shù)的估計精度有明顯的提升，但其精度受限于n的取值，只能得到1/n精度估計結(jié)果。而本文算法在1/n精度估計的基礎(chǔ)上，利用擬合修正估計值，實現(xiàn)了任意精度估計結(jié)果。

實驗中，通過計算可得，F(xiàn)M_sinc，F(xiàn)ast_FM和本文算法的平均運(yùn)行時間分別為0.763 7 s，0.790 1 s和0.805 3 s，算法運(yùn)行時間分別增加約0.041 6 s和0.015 2 s。說明本文算法能在增加較少運(yùn)行時間的情況下獲得更高的配準(zhǔn)精度。

仿真圖像在平移固定，旋轉(zhuǎn)縮放(15°，1.2)下，采用3種算法配準(zhǔn)后的圖像與參考圖像之間的聯(lián)合概率直方圖如圖4所示。

圖4 旋轉(zhuǎn)縮放(15°，1.2)時，3種算法配準(zhǔn)結(jié)果的聯(lián)合直方圖

從圖4可以看出：圖4(a)中有大量的點分散在對角線周圍，相比圖4(a)，圖4(b)中有部分點更靠近對角線，圖4(c)中點的分布最靠近直角線。聯(lián)合直方圖中的點分布越靠近對角線，誤差越小，圖4結(jié)果表明本文算法的配準(zhǔn)誤差最小。

3.2 真實圖像實驗

圖5 真實圖像

選取某地區(qū)2001年拍攝的SPOT4圖像和2004年拍攝的Landsat TM圖像作為真實圖像進(jìn)行實驗。分別從中心截取大小為330×330的子圖作為參考圖和待配準(zhǔn)圖，如圖5所示。

由于真實遙感圖像間的旋轉(zhuǎn)、縮放和平移量未知，無法對結(jié)果客觀評價，所以根據(jù)配準(zhǔn)后的拼接效果進(jìn)行主觀評價。圖6給出3種算法配準(zhǔn)后的拼接結(jié)果圖，其中每幅圖右邊的小圖依次為左圖中從左到右的圓圈中拼接區(qū)域的局部放大圖。從圖6可以看出，3種算法都得到了很好的配準(zhǔn)結(jié)果，但仔細(xì)分析各圖的第一個局部放大區(qū)域可以看出，F(xiàn)M_sinc求得的配準(zhǔn)后圖像在拼接處有完全錯開的斷層，F(xiàn)ast_FM算法求得的圖像拼接處斷層錯開近1/4，本文算法比Fast_FM的拼接效果稍微好一些，斷層處錯開了近1/2。由此可見，本文算法的配準(zhǔn)效果最佳。

圖6 3種算法配準(zhǔn)拼接結(jié)果

4 結(jié)束語

針對亞像素級的旋轉(zhuǎn)、縮放運(yùn)動的高精度估計問題，本文提出一種基于曲面擬合修正矩陣乘法的傅里葉梅林變換算法，計算結(jié)果表明本文算法的計算精度誤差為0.15像素級別，是一種計算精度高的配準(zhǔn)方法，具有廣泛的應(yīng)用前景。如何進(jìn)一步提高計算精度，增強(qiáng)其魯棒性，減少計算時間是下一步研究的重點。