焦繼超，趙保軍，周剛

(1.北京理工大學(xué) 信息與電子學(xué)院，北京100081;2.中國(guó)科學(xué)院 長(zhǎng)春光學(xué)精密機(jī)械與物理研究所，吉林 長(zhǎng)春130033)

0 引言

為了檢測(cè)跟蹤空間中的暗弱碎片，需要口徑大而焦距短的望遠(yuǎn)鏡。但是，隨著望遠(yuǎn)鏡口徑的擴(kuò)大，望遠(yuǎn)鏡的焦距也要變長(zhǎng)，于是望遠(yuǎn)鏡的CCD 視場(chǎng)就會(huì)變小［1］，這就降低了地基天文望遠(yuǎn)鏡探測(cè)空間碎片的能力。同時(shí)，研制大口徑的天文望遠(yuǎn)鏡也存在技術(shù)、研制周期等方面的困難。此時(shí)，采用并行排列的望遠(yuǎn)鏡陣列能夠解決上述問(wèn)題［2］。由于采用陣列望遠(yuǎn)鏡在同一時(shí)刻拍攝的圖像并不完全重合，需要對(duì)不同CCD 相機(jī)拍攝的圖像進(jìn)行配準(zhǔn)，才能實(shí)現(xiàn)圖像的精確疊加，提高圖像信噪比，以檢測(cè)出暗弱的空間碎片。

基于區(qū)域的方法是從參考圖像中選取一定大小的區(qū)域在待配準(zhǔn)圖像中搜索與其相關(guān)性最大的圖像塊，把2 個(gè)區(qū)域圖像的中心作為一對(duì)特征點(diǎn)，然后利用這些特征點(diǎn)求解變換方程［3］。利用灰度值作為相似度測(cè)量的算法結(jié)構(gòu)簡(jiǎn)單［4］，易于實(shí)現(xiàn)，但是運(yùn)算量大，容易受光照、噪聲等因素的干擾。采用互信息配準(zhǔn)的方法通過(guò)計(jì)算圖像概率密度的相關(guān)性實(shí)現(xiàn)圖像的精確配準(zhǔn)，多用于醫(yī)學(xué)圖像處理中。但是該算法需要估計(jì)圖像的概率密度，因此對(duì)于概率分布比較復(fù)雜的圖像具有較大的配準(zhǔn)誤差。同時(shí)，空間圖像中的星體基本都呈斑狀，結(jié)構(gòu)特征較少，采用基于結(jié)構(gòu)特征的配準(zhǔn)方法存在較大困難。

針對(duì)一般算法難以實(shí)現(xiàn)大分辨率、低信噪比空間圖像實(shí)時(shí)和精確配準(zhǔn)的問(wèn)題，本文采用了一種基于改進(jìn)的傅里葉-梅林變換的配準(zhǔn)算法。實(shí)驗(yàn)結(jié)果表明本算法能夠?qū)崿F(xiàn)空間圖像快速配準(zhǔn)的前提下，獲得亞像素的配準(zhǔn)精度。

1 基于FMT 的圖像配準(zhǔn)基本原理

基于FMT 的配準(zhǔn)算法是從基于FFT 的配準(zhǔn)算法發(fā)展起來(lái)的，它利用圖像頻譜的相位相關(guān)性實(shí)現(xiàn)存在平移、旋轉(zhuǎn)和縮放的2 幅圖像間的配準(zhǔn)［5］。

對(duì)于存在平移、旋轉(zhuǎn)和縮放的圖像f1(x，y)和f2(x，y)，兩者在空域中的位置關(guān)系如下所示

其中，f1為參考圖像，f2為待配準(zhǔn)圖像，為了便于敘述，圖像在水平和垂直方向采用相同的縮放系數(shù)a，θ0為圖像的旋轉(zhuǎn)系數(shù)，Δx 和Δy 分別為水平方向和垂直方向上的平移量。

對(duì)式(1)進(jìn)行傅里葉變換［6］，得到如下等式

其中:F1(ξ，η)為f1的頻譜，F(xiàn)2(ξ，η)為f2的頻譜。由上式可知，圖像在傅里葉變換前后的旋轉(zhuǎn)角度不變，縮放系數(shù)變?yōu)樵档牡箶?shù)。不考慮圖像間的平移運(yùn)動(dòng)時(shí)，F(xiàn)1和F2的幅度值間的關(guān)系如下

由上式可知，方程具有三角函數(shù)的形式，因此可以將式(3)中的變量進(jìn)行如下形式的代換

其中，rp(θ，ρ)和sp(θ，ρ)分別為F1和F2在極坐標(biāo)中的頻譜。需要注意的是，頻譜的幅度值對(duì)應(yīng)的旋轉(zhuǎn)角度是周期性的，其周期為π，則rp和sp有如下性質(zhì)［7］

其中，n=…-2，-1，0，1，2，….則式(3)中的坐標(biāo)變量進(jìn)行極坐標(biāo)轉(zhuǎn)換，其形式如下所示

其中，θ 為極坐標(biāo)系下的角度參數(shù)，θ0為圖像間的旋轉(zhuǎn)角度。把式(8)和式(9)帶入(3)式可得

則上式可以簡(jiǎn)化為

由上式可知，在極坐標(biāo)系下，圖像旋轉(zhuǎn)角度θ 沿極點(diǎn)進(jìn)行旋轉(zhuǎn)，縮放系數(shù)在極軸方向上成比例變化。因此，對(duì)變量進(jìn)行對(duì)數(shù)變換

其中:λ 為極坐標(biāo)系下的極徑ρ 的對(duì)數(shù)，b 為直角坐標(biāo)系下縮放系數(shù)a 的對(duì)數(shù)。同時(shí)我們定義:rpl(θ，λ)=rp(θ，ρ)，spl(θ，λ)=sp(θ，ρ)，則式(10)可以變?yōu)槿缦滦问?/p>

由式(13)可知，在對(duì)數(shù)極坐標(biāo)系下，不但圖像旋轉(zhuǎn)角度θ 沿極點(diǎn)進(jìn)行旋轉(zhuǎn)，而且縮放系數(shù)a 也沿極軸進(jìn)行平移。這2 個(gè)性質(zhì)被稱為角度不變性和距離不變性?；谙辔幌嚓P(guān)的方法獲得旋轉(zhuǎn)系數(shù)和縮放系數(shù)，則式(13)兩端進(jìn)行傅里葉變換，

把上式帶入交叉能量譜公式中［8］，

當(dāng)?shù)仁阶筮呌凶畲笾禃r(shí)，求得的θ0和b 即為在對(duì)數(shù)極坐標(biāo)系中對(duì)應(yīng)的旋轉(zhuǎn)系數(shù)和縮放系數(shù)。因此可以根據(jù)式(16)求得在直角坐標(biāo)系下的縮放系數(shù)

2 新的基于FMT 的空間圖像配準(zhǔn)算法的原理及實(shí)現(xiàn)

2.1 一種新的旋轉(zhuǎn)系數(shù)判斷方法

空間圖像的分辨率一般在1 024 ×1 024 像素或2 048 ×2 048 像素甚至更高，即使采用快速傅里葉變換也要進(jìn)行大量的計(jì)算，影響了空間圖像配準(zhǔn)的實(shí)時(shí)性。同時(shí)，從第二部分對(duì)基于FMT 配準(zhǔn)算法原理的論述來(lái)看，通過(guò)提高圖像傅里葉變換的速度來(lái)保持配準(zhǔn)算法的實(shí)時(shí)性比較困難，因?yàn)槟壳斑€沒(méi)有比FFT 變換更為簡(jiǎn)單快速的傅里葉變換算法。一方面，拍攝得到的空間圖像間存在的平移運(yùn)動(dòng)和縮放對(duì)空間圖像的相關(guān)性影響相對(duì)較小;另一方面，在頻域中，由于旋轉(zhuǎn)系數(shù)存在周期性，為了獲得正確的系數(shù)，需要利用獲得的縮放系數(shù)和多個(gè)旋轉(zhuǎn)系數(shù)對(duì)待配準(zhǔn)圖像進(jìn)行轉(zhuǎn)換并計(jì)算其與參考圖像的相關(guān)系數(shù)，這一過(guò)程需要大量的計(jì)算。

為了提高空間圖像配準(zhǔn)的速度，保證其實(shí)時(shí)性，本文提出了一種新的旋轉(zhuǎn)系數(shù)判斷方法，即通過(guò)計(jì)算兩幅原始配準(zhǔn)圖像間的相關(guān)值T，并與閾值βth進(jìn)行比較，以獲得正確的旋轉(zhuǎn)系數(shù)，其具體的判斷步驟如下:

1)根據(jù)相位相關(guān)性，得到在對(duì)數(shù)極坐標(biāo)系下配準(zhǔn)圖像間的旋轉(zhuǎn)系數(shù)θ1和θ2且θ1＞θ2;

2)計(jì)算參考空間圖像和待配準(zhǔn)空間圖像間的相關(guān)值T，公式如下所示［9］，

其中，M 為圖像列向量的個(gè)數(shù)，N 為圖像行向量的個(gè)數(shù)，fs(x，y)為待配準(zhǔn)圖像，fr(x，y)為參考圖像;

3)比較T 與閾值βth間的大小，如果T ＞βth，則說(shuō)明兩幅圖像間的相關(guān)性較大，選取θ2為旋轉(zhuǎn)系數(shù);如果T＜βth，則說(shuō)明兩幅圖像間的差別比較大，選取θ1為旋轉(zhuǎn)系數(shù)，其中βth為常數(shù)，由大量的實(shí)驗(yàn)得到，本文試驗(yàn)中βth取值范圍是3 500～3 900.

從以上的步驟來(lái)看，本文提出的方法只進(jìn)行2次相關(guān)運(yùn)算，而原旋轉(zhuǎn)角的判斷方法至少要進(jìn)行4次乘法運(yùn)算(包括對(duì)圖像放大σ-1倍時(shí)需要插值的2 次運(yùn)算及計(jì)算交叉能量譜的值)和計(jì)算交叉能量譜時(shí)對(duì)圖像進(jìn)行的2 次FFT 變換。以1 024 ×1 024像素的空間圖像為例，2 種方法中主要算法的運(yùn)算量如表1所示。

表1 2 種旋轉(zhuǎn)系數(shù)判斷方法運(yùn)算量的比較Tab.1 The comparison of computation between two kinds of rotation coefficient determination methods

表中，表中的除法是計(jì)算交叉能量譜時(shí)用到的。由表1可以看出，本文提出的判斷方法之所以能夠減少運(yùn)算量，增加配準(zhǔn)算法的運(yùn)算速度，就是因?yàn)榇蟠鬁p少了計(jì)算結(jié)構(gòu)復(fù)雜的除法和復(fù)數(shù)乘法的次數(shù)。

2.2 空間圖像配準(zhǔn)算法中的插值

本文提出的空間圖像配準(zhǔn)算法需要在不同階段進(jìn)行多次插值，采用相同的插值方法會(huì)影響整個(gè)配準(zhǔn)算法的運(yùn)行速度，因此本文在配準(zhǔn)的不同階段采用不同的插值方法，在不降低配準(zhǔn)速度的前提下實(shí)現(xiàn)較高的配準(zhǔn)精度。

在基于FMT 空間圖像配準(zhǔn)算法中，為了計(jì)算旋轉(zhuǎn)系數(shù)和縮放系數(shù)，需要進(jìn)行極坐標(biāo)和對(duì)數(shù)坐標(biāo)的轉(zhuǎn)換，極坐標(biāo)系和對(duì)數(shù)極坐標(biāo)系的坐標(biāo)點(diǎn)都是非均勻分布的，而直角坐標(biāo)系坐標(biāo)點(diǎn)是均勻分布的，直角坐標(biāo)系中的所有點(diǎn)在轉(zhuǎn)換后不一定都對(duì)應(yīng)其它兩種坐標(biāo)系中的整數(shù)點(diǎn)，因此需要進(jìn)行插值。因?yàn)樵谵D(zhuǎn)換的過(guò)程中需要進(jìn)行2 次插值，為了保證空間圖像配準(zhǔn)的實(shí)時(shí)性，需要采用易于實(shí)現(xiàn)的插值方法，本文采用的是雙線性插值法。

一般情況下，計(jì)算得到的轉(zhuǎn)換參數(shù)是非整數(shù)，當(dāng)利用轉(zhuǎn)換方程對(duì)待配準(zhǔn)圖像進(jìn)行轉(zhuǎn)換時(shí)，輸出像素通常被映射到新坐標(biāo)系下的非整數(shù)位置，因此為了決定與該位置相對(duì)應(yīng)的灰度級(jí)，也要進(jìn)行插值。這個(gè)階段的插值對(duì)最后的配準(zhǔn)精度有著較大影響，所以采用3 次樣條插值。常用的插值方法主要有最鄰近插值、雙線性插值、樣條插值。最鄰近插值和雙線性插值簡(jiǎn)單易于實(shí)現(xiàn)，但是人工處理的痕跡比較明顯，3 次樣條插值能夠較好的消除鋸齒現(xiàn)象，插值精度比較高［10］，3 次樣條插值的具體步驟如下，

1)輸入m 個(gè)插值結(jié)點(diǎn)，α =x1＜x2＜…＜xm=β，對(duì)應(yīng)函數(shù)值為y1，y2，…，ym，邊界條件y'1，y'2，…，y'm，待求插值點(diǎn)x0;

2)計(jì)算gi=xi+1-xi，其中i=1，2，…，m-1;

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

5)在滿足邊界條件的前提下，計(jì)算如下方程

其中，φn和ωn為插值基函數(shù)，

6)輸出各區(qū)間的3 次樣條插值函數(shù)si(x).

2.3 空間圖像配準(zhǔn)的流程圖

綜合前面的論述，圖1給出了本文提出的算法對(duì)參考空間圖像s(x，y)和待配準(zhǔn)空間圖像r(x，y)進(jìn)行配準(zhǔn)的流程圖，其中θ1＞θ2，如下所示，

圖1 空間圖像配準(zhǔn)流程圖Fig.1 The space image registration flowchart

3 實(shí)驗(yàn)結(jié)果與分析結(jié)論

在所研制的DSP+FPGA 硬件平臺(tái)上通過(guò)了實(shí)驗(yàn)驗(yàn)證，結(jié)構(gòu)框圖如圖2所示，其中，實(shí)地拍攝的空間圖像在 FPGA 中進(jìn)行預(yù)處理，在 DSP(TMS320C6455)中進(jìn)行基于FMT 的圖像配準(zhǔn)，DPRAM 存儲(chǔ)原始空間圖像，SSRAM 緩存圖像數(shù)據(jù)，處理結(jié)果發(fā)送給主控計(jì)算機(jī)。

圖2 DSP+FPGA 結(jié)構(gòu)框圖Fig.2 Block Diagram of DSP+FPGA

為了驗(yàn)證本算法的性能，本文采用實(shí)地拍攝且分辨率為1 024 ×1 024 像素的空間圖像為參考圖像r(x，y)，對(duì)r(x，y)進(jìn)行向下和向右平移(Δx =2，Δy=1)，并作3°的逆時(shí)針旋轉(zhuǎn)得到待配準(zhǔn)圖像s(x，y)，如圖3所示，2 幅圖像的噪聲比較嚴(yán)重，其中標(biāo)記出的是待配準(zhǔn)圖像的獨(dú)有星體。

為了證明本算法在空間圖像配準(zhǔn)中的性能，不但采用均方根誤差(RMSE)準(zhǔn)則對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行評(píng)價(jià)，而且還同基于區(qū)域配準(zhǔn)和基于結(jié)構(gòu)特征的配準(zhǔn)算法進(jìn)行比較。

圖3 配準(zhǔn)空間圖像Fig.3 The registration space images

3.1 本文配準(zhǔn)算法性能分析

由圖3可以發(fā)現(xiàn)，空間圖像的信噪比很低，有著很多明顯的噪點(diǎn)，但是通過(guò)高通濾波可增強(qiáng)星體邊緣，抑制部分噪聲，圖4給出了圖3經(jīng)過(guò)傅里葉變換和濾波后的三維頻譜圖。

圖4 空間圖像直角坐標(biāo)系頻譜圖Fig.4 The space image spectrum in cartesian coordinate

從圖4可以看出，空間圖像頻譜中某一個(gè)頻率對(duì)應(yīng)一個(gè)全局的峰值，而噪聲對(duì)應(yīng)的幅度值要小于這個(gè)值，因此可利用2 個(gè)頻譜幅度值峰值的相關(guān)性實(shí)現(xiàn)圖像的配準(zhǔn)。

為了獲得配準(zhǔn)轉(zhuǎn)換方程的旋轉(zhuǎn)系數(shù)和縮放系數(shù)，需要將直角坐標(biāo)系下的頻譜轉(zhuǎn)換到對(duì)數(shù)極坐標(biāo)系下，圖5給出了對(duì)數(shù)極坐標(biāo)系下的空間圖像頻譜圖。

由圖5可知，對(duì)數(shù)極坐標(biāo)系下的空間圖像頻譜和直角坐標(biāo)系下的頻譜類似，某一坐標(biāo)值對(duì)應(yīng)一個(gè)全局最大幅度值。

利用相位相關(guān)性，分別在對(duì)數(shù)極坐標(biāo)系和直角坐標(biāo)系中求出旋轉(zhuǎn)角度θ =3.02°、縮放系數(shù)a =1.01 和平移系數(shù)(Δx =2.00，Δy =1.00)，圖6給出了配準(zhǔn)后的圖像s'(x，y)以及和r(x，y)的疊加圖像c(x，y).

從疊加圖像c(x，y)中可以看出，不論是面積很小的星體還是拖尾都能夠完全重合，而沒(méi)有明顯的疊加痕跡，其中圓形標(biāo)記的和圖3(b)中標(biāo)記出的是同一星體。

同時(shí)本文還采用均方根誤差(RMSE)準(zhǔn)則來(lái)評(píng)價(jià)算法性能，其定義如下［11］

圖5 空間圖像對(duì)數(shù)極坐標(biāo)系頻譜圖Fig.5 The space image spectrum in log-polar coordinate

其中，N 為空間圖像提取的特征星體個(gè)數(shù)，(x，y)為參考圖像中特征星體的質(zhì)心位置，(x'，y')為配準(zhǔn)后圖像中特征星體的質(zhì)心位置。

利用式(19)可以求得空間圖像的RMSE，配準(zhǔn)前，RMSE=1.202，配準(zhǔn)后，RMSE=0.518，即配準(zhǔn)前后RMSE 減少值ΔRMSE =0.684.從主觀觀測(cè)和客觀指標(biāo)方面都說(shuō)明了本文配準(zhǔn)算法對(duì)空間圖像進(jìn)行了精確配準(zhǔn)。

3.2 本文算法與其它配準(zhǔn)算法的性能比

為了進(jìn)一步突出本算法的性能，本次試驗(yàn)又利用基于區(qū)域配準(zhǔn)算法和基于結(jié)構(gòu)特征提取配準(zhǔn)算法進(jìn)行空間圖像的配準(zhǔn)，表2在運(yùn)算時(shí)間和RMSE 兩個(gè)指標(biāo)上對(duì)3 種配準(zhǔn)算法進(jìn)行了比較。

由表2可知，在運(yùn)算速度上，本文提出的算法采用了新的旋轉(zhuǎn)角度判斷方法和多種插值方式，因此分別比基于區(qū)域和基于結(jié)構(gòu)特征的配準(zhǔn)算法分別快4.155 s 和3.476 s;在RMSE 上，基于FMT 的配準(zhǔn)算法比其他兩種配準(zhǔn)算法分別減少0.436 s 和0.188 s.

圖6 配準(zhǔn)后和疊加空間圖像Fig.6 The space image after registration and superposable image

表2 3 種配準(zhǔn)算法的性能比較Tab.2 The comparison of the three registration methods

4 結(jié)論

本文根據(jù)空間圖像配準(zhǔn)實(shí)時(shí)性和精確性的要求，提出了一種新的基于FMT 的配準(zhǔn)算法。首先，提出一種新的基于原始空間圖像相關(guān)性的旋轉(zhuǎn)系數(shù)判斷方法。然后，為了在不損失配準(zhǔn)精度前提下降低算法復(fù)雜度，提出在坐標(biāo)系轉(zhuǎn)換和待配準(zhǔn)圖像變換階段采用不同插值方法，并給出了三次樣條插值的步驟。最后，給出了本文算法實(shí)現(xiàn)空間圖像配準(zhǔn)的流程圖。通過(guò)硬件平臺(tái)的實(shí)驗(yàn)驗(yàn)證，結(jié)果證明:基于FMT 的配準(zhǔn)方法能夠基本實(shí)現(xiàn)空間圖像配準(zhǔn)的實(shí)時(shí)性要求;配準(zhǔn)精度達(dá)到0.5 像素;算法的RMSE減少0.684.基本滿足了空間圖像配準(zhǔn)對(duì)實(shí)時(shí)性和精確性的要求。

References)

［1］王鳴浩，陳濤，王建立，等.捆綁式望遠(yuǎn)鏡圖像信噪比測(cè)量及分析［J］.光學(xué)精密工程，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］劉衛(wèi)光，崔江濤，周利華.插值和相位相關(guān)的圖像亞像素配準(zhǔn)方法［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，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 算法的干涉圖像配準(zhǔn)［J］.光學(xué)精密工程，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.