李喜艷，紀(jì)東升，吳崇正

(1．鄭州成功財(cái)經(jīng)學(xué)院 信息工程系，河南 鄭州451200;2．蘭州大學(xué) 信息科學(xué)與工程學(xué)院，甘肅 蘭州730030)

0 引言

圖像拼接是通過對(duì)多個(gè)具有互相重疊部分的圖像進(jìn)行配準(zhǔn)，拼接成一個(gè)分辨率更高，視野更寬的圖像．圖像配準(zhǔn)是圖像拼接的重要步驟，基于特征點(diǎn)匹配的方法須先從圖像上提取不變性特征，如特征角點(diǎn)、矩等，再進(jìn)行精確匹配．在圖像的各種不變性特征中，角點(diǎn)具有旋轉(zhuǎn)不變性和幾乎不受光照條件影響的優(yōu)點(diǎn)，因此，基于特征點(diǎn)匹配的圖像拼接在視頻檢索、景物匹配以及遙感圖像融合等領(lǐng)域都具有十分廣泛的應(yīng)用． 為有效解決圖像在旋轉(zhuǎn)、平移變換及光度差異等情況下存在的拼接問題，提出了一種基于改進(jìn)的SURF 檢測(cè)［1］和偽Zernike 矩的圖像拼接算法(SURF Feature Points Pseudo-Zernike Moments Algorithm，SPZMA)．該算法利用形態(tài)學(xué)極大值改進(jìn)的SURF 算子提取圖像特征點(diǎn);然后計(jì)算中心鄰域窗口的偽Zernike 矩，通過比較偽Zernike 矩的Bray-Curtis距離得到初始匹配點(diǎn)對(duì);接著運(yùn)用RANSAC 算法［2］剔除偽特征點(diǎn)對(duì)，建立圖像對(duì)之間的變換模型，實(shí)現(xiàn)圖像配準(zhǔn);最后，通過利用漸入漸出分段加權(quán)函數(shù)消除拼接縫，使拼接圖像實(shí)現(xiàn)平穩(wěn)過渡，SPZMA 算法流程結(jié)構(gòu)如圖1 所示．

圖1 SPZMA 算法流程圖Fig．1 SPZMA algorithm flowchart

1 圖像特征點(diǎn)檢測(cè)和描述

1．1 SURF 特征點(diǎn)檢測(cè)算法

(1)積分圖像中任意一點(diǎn)(i，j)的值ii(i，j)為原圖像左上角到任意點(diǎn)(i，j)相應(yīng)的對(duì)角線區(qū)域灰度值的總和，即

式中:p(m，n)表示原圖像中點(diǎn)(m，n)的灰度值．

(2)給定圖像I 中一個(gè)點(diǎn)x(x，y)，在點(diǎn)x 處，尺度為σ 的Hessian 矩陣H(x，σ)定義如下:

式中:Lxx(x，σ)是高斯二階微分(σ)在點(diǎn)x處與圖像I 的卷積，Lxy(x，σ)和Lyy(x，σ)具有同樣的含義．

為了將模板與圖像的卷積轉(zhuǎn)化成濾波［2］運(yùn)算，需要對(duì)高斯二階微分模板進(jìn)行簡(jiǎn)化，使簡(jiǎn)化后的模板只由幾個(gè)矩形區(qū)域組成，矩形區(qū)域內(nèi)填充同一值，用Dxx、Dyy和Dxy表示模板與圖像卷積的結(jié)果，這樣便可將Hessian 矩陣［3］的行列式作如下簡(jiǎn)化，其中w 為權(quán)重．

Det(Happrox)=DxxDyy-(wDxy)2． (3)

(3)將二維圖像I (x，y)與高斯核卷積得到不同尺度的空間:L(x，y，σ ) = I (x，y)·G(x，y，σ )，其中，G(x，y，σ )為尺度可變的高斯函數(shù)，再通過建立高斯差分尺度空間尋找圖像的局部極值．在圖像高斯差分尺度空間內(nèi)當(dāng)前尺度和其相鄰兩個(gè)尺度3* 3 的區(qū)域內(nèi)，標(biāo)記的像素點(diǎn)和周圍像素點(diǎn)比較，如果標(biāo)記的像素灰度值大于或者小于其他像素，那么這個(gè)標(biāo)記像素點(diǎn)就是特征點(diǎn)．

Step4 對(duì)上一步的小數(shù)組做2(w -1)次比較．這里對(duì)算法進(jìn)行如下改進(jìn):

計(jì)算Rk和Qk平均需要進(jìn)行max(w -q，q -1)+(w-1)+1 =1．5w -(w mod2)/2 次計(jì)算和比較，且小于2(w-1)．

Step5 窗口最大值的函數(shù)可表示為

1．2 利用局部極值搜索窗口算法的改進(jìn)

根據(jù)Hessian 矩陣求出尺度圖像在(x，y，σ)的極值后，首先在極值點(diǎn)的3 ×3 ×3 的立體領(lǐng)域內(nèi)進(jìn)行非極大值抑制［4］． 為了能夠?qū)蜻x特征點(diǎn)進(jìn)行亞像素定位，可以在尺度空間和圖像空間中使用二次擬合函數(shù)進(jìn)行插值:

筆者基于形態(tài)學(xué)濾波中的局部極大值抑制的思想，采用改進(jìn)的HGW(Herk-Gil-Werman)算法，將局部求極大值抑制的過程轉(zhuǎn)化為在二維窗口內(nèi)對(duì)局部極值的求取，用以加快算法的速度．

Step1 先將像素矩陣為X 的圖像擴(kuò)展為Y=［Y A］，之后將擴(kuò)展矩陣均分為寬度為w 的多個(gè)子矩陣，A 是擴(kuò)展的零矩陣．

Step2 將二維窗口劃分為2 ×1． 極值二維窗口用下式定義為

Step3 確定擴(kuò)展矩陣Y 某行數(shù)組(x0，…，xn-1)的局部最大值． 用公式y(tǒng)i=xi+j(i =0，…，n-w)表示求最大值，且將x0，…，xn-1分成多個(gè)小數(shù)組xw-1，x2w-1，x3w-1，…．對(duì)每個(gè)小數(shù)組，定義如下序列元素Rk和QK(k=0，…，w-1 )．

應(yīng)用二叉樹搜索定理可知，每個(gè)元素需要搜索「lb(w -1)?/w 次，進(jìn)而可以在擴(kuò)展矩陣Y 內(nèi)確定最大值的新矩陣Z，故數(shù)組內(nèi)每個(gè)元素平均比較次數(shù)為

1．5 +「lb(w-1)?/w-(w mod2)/(2w)． (9)

Step6 重復(fù)進(jìn)行上述算法的Step1 至Step5，即可得到二維窗口局部極值的分布矩陣Z～，即Z～表示在某元素半徑范圍內(nèi)的最大值的分布情況．

Step7 去掉分布矩陣Z～中的擴(kuò)充矩陣，且與X 矩陣對(duì)應(yīng)元素作減法，若有0 值，便為w ×w 窗口的局部最大值．

由于在垂直與水平方向上均進(jìn)行了一次局部極值的求解，改進(jìn)后的總的元素比較次數(shù)就為

3 +2「lb(w-1)?/w-(w mod2)/w． (10)

圖2 為改進(jìn)后的HGW 算法、原HGW 算法與普通算法實(shí)驗(yàn)結(jié)果對(duì)比．可以看到，改進(jìn)后的算法單個(gè)元素的比較次數(shù)將趨向于3，比其他兩種算法都要少．當(dāng)窗口越來越大時(shí)，改進(jìn)算法的性能優(yōu)越性就越來越明顯．

圖2 3 種算法的實(shí)驗(yàn)效果比較Fig．2 Experimental results compare the three algorithm

實(shí)驗(yàn)表明，改進(jìn)后的HGW 算法角點(diǎn)檢測(cè)定位準(zhǔn)確，自適應(yīng)性好，角點(diǎn)提取速度更快，從而得到比原算法更理想的角點(diǎn)檢測(cè)效果．圖3 ～6 為本文算法對(duì)角點(diǎn)的偏移、lena 圖像、旋轉(zhuǎn)45°以及加入噪聲的檢測(cè)效果．

圖3 表明，改進(jìn)算法定位關(guān)鍵點(diǎn)角點(diǎn)更加準(zhǔn)確．圖4 和圖5 可看出，本文算法提取的角點(diǎn)均勻，對(duì)圖像旋轉(zhuǎn)的適應(yīng)性好．圖6 中加入了乘性噪聲(σ=0．02)，本文算法和原算法相比較，提取角點(diǎn)均勻合理，具有較強(qiáng)的抗噪能力．表1 為圖像角點(diǎn)提取個(gè)數(shù)及平均運(yùn)行時(shí)間的統(tǒng)計(jì)，運(yùn)用快速窗口搜索算法的改進(jìn)算法在運(yùn)行速度上有了明顯提高，雖然提取出的角點(diǎn)數(shù)量比原算法減少，但角點(diǎn)分布更加均勻，尤其對(duì)關(guān)鍵點(diǎn)角點(diǎn)檢測(cè)沒有漏檢．

圖3 本文改進(jìn)算法與原算法對(duì)克服角點(diǎn)偏移結(jié)果比較Fig．3 The improved algorithm with the original algorithm compare the result to overcome offset corner

圖4 本文改進(jìn)算法與原算法對(duì)Lena 圖像實(shí)驗(yàn)結(jié)果比較Fig．4 The experimental results on the Lena image rotation 45 degrees compare

圖5 對(duì)Lena 圖像旋轉(zhuǎn)45 度的實(shí)驗(yàn)結(jié)果比較Fig．5 The experimental results on the Lena image rotation 45 degrees compare

圖6 加入乘性噪聲(σ=0．02)的實(shí)驗(yàn)結(jié)果比較Fig．6 The join multiplicative noise(σ=0．02)of the experimental results

表1 檢測(cè)Lena 圖像角點(diǎn)個(gè)數(shù)及平均運(yùn)行時(shí)間統(tǒng)計(jì)Tab．1 Lena image corner detection of the number of points and the average run time statistics

1．3 特征點(diǎn)描述

(1)特征點(diǎn)提取后，需要找出兩幅圖像特征點(diǎn)之間的對(duì)應(yīng)關(guān)系，即特征匹配． Zernike 矩是一組正交矩，具有旋轉(zhuǎn)不變性的特性，且Zernike 矩可以構(gòu)造任意高階矩，所以識(shí)別效果優(yōu)于其他方法．筆者采用特征點(diǎn)鄰域的不變矩作為匹配的特征，根據(jù)偽Zernike 矩的不變性，使它具有更多的矩?cái)?shù)量及更好的抗噪聲性能．

(2)偽Zernike 矩的選擇．?dāng)?shù)字圖像的像素離散性質(zhì)使得偽Zernike 矩的計(jì)算產(chǎn)生誤差，不同階矩的計(jì)算精確度不同［5］，因此，必須對(duì)矩進(jìn)行優(yōu)化選擇．選擇矩時(shí)主要考慮兩點(diǎn):1)當(dāng)矩的階數(shù)高于某一個(gè)值Nmax時(shí)，計(jì)算量大且計(jì)算不再準(zhǔn)確，實(shí)驗(yàn)中取10 階;2)重復(fù)度為m =4i(i =0，1，2，…．)的矩計(jì)算是不準(zhǔn)確的，應(yīng)當(dāng)去除． 而由于矩的共扼對(duì)稱性，剩下的矩只有一半是獨(dú)立的．記選取的矩集合為S，即S={Anm，n≤Nmax，m≥0，m≠4i}．

2 圖像匹配

2．1 圖像的快速匹配

為了達(dá)到好的特征點(diǎn)匹配效果，采用“粗—精”的方法以提高特征點(diǎn)匹配的可靠性［6］． ①用Bray-Curtis 相似性度量初步選取偽Zernike 矩描述的相似匹配特征點(diǎn)對(duì);②采用RANSAC 算法進(jìn)一步篩選偽匹配點(diǎn)對(duì)．

歐式距離度量通常被用于計(jì)算圖像之間的相似度，然而，其不能對(duì)特征值進(jìn)行歸一化，會(huì)有一些高量級(jí)的特征點(diǎn)控制著低量級(jí)的特征點(diǎn)． 筆者運(yùn)用Bray-Curtis 相似性度量，規(guī)范其特征值并通過區(qū)分不同的相對(duì)特征值和絕對(duì)值，以計(jì)算相似特征點(diǎn)對(duì)的匹配度． Bray-Curtis 的相似性度量公式(11)提出了基于特征點(diǎn)的描述:

式中:fi(T)、fi(D)表示特征向量的測(cè)試和數(shù)據(jù)圖像;H 是基于角點(diǎn)特征的數(shù)量．

Bray-Curtis 的相似性度量用于偽Zernike 矩:

d(T，D)=wcdc(T，D)+wrdr(T，D)． (13)

式中:wc和wr代表權(quán)重的兩圖像的相似性度量．計(jì)算相似度和Bray-Curtis 相似性度量?jī)蓚€(gè)歐氏空間，比較它們的性能．這樣可以計(jì)算出原圖像中每一個(gè)特征點(diǎn)在目標(biāo)圖像中的匹配點(diǎn)，當(dāng)然這樣的匹配對(duì)中會(huì)包含部分偽匹配點(diǎn)對(duì)．

2．2 精確匹配

RANSAC 算法根據(jù)一組包含異常數(shù)據(jù)的樣本數(shù)據(jù)集，計(jì)算出數(shù)據(jù)的數(shù)學(xué)模型參數(shù)，并得到有效的樣本數(shù)據(jù)．利用齊次坐標(biāo)，兩幅圖像之間的投影變換模型可以用矩陣的形式來描述［7］:

式中:參數(shù)m11，m12，m21，m22表示尺度縮放和旋轉(zhuǎn);m13為水平方向位移;m23為垂直方向位移．

以相鄰兩幅圖像為例，設(shè)匹配特征點(diǎn)對(duì)數(shù)目為N，匹配特征點(diǎn)集合分別記為P(1，N)，P(2，N)．其中，P(1，N)為基準(zhǔn)圖像的特征點(diǎn)集合，P(2，N)為待匹配圖像的特征點(diǎn)集合，算法步驟如下:

Step1 從初始N 對(duì)匹配特征點(diǎn)中隨機(jī)選取3 對(duì)匹配特征點(diǎn);

Step2 由選取的3 對(duì)匹配特征點(diǎn)計(jì)算出基準(zhǔn)圖像和待匹配圖像間的相似變換矩陣M12，利用變換矩陣M12 將待匹配圖像的特征點(diǎn)集合P(2，N)中剩余的N-2 個(gè)特征點(diǎn)P(2，N -2)變換到基準(zhǔn)圖像坐標(biāo)系下，并記為P'(2，N-2);

Step3 計(jì)算變換后的特征點(diǎn)P'(2，N-2)與特征點(diǎn)P(1，N-2)之間的坐標(biāo)誤差;

Step4 從N 對(duì)匹配特征點(diǎn)對(duì)中找出坐標(biāo)誤差在一定誤差閾值內(nèi)特征點(diǎn)對(duì)個(gè)數(shù)，記為i;

Step5 反復(fù)迭代Step1 ～Step4 步，直至找到i 最大的集合為最大內(nèi)點(diǎn)集即為內(nèi)點(diǎn)，其余N -i為誤匹配點(diǎn)即為外點(diǎn)．用最小二乘法來減少誤差，去除了誤匹配的影響，最終得到空間變換矩陣M．

3 圖像融合

采用具有雙線插值思想的漸入漸出方法［8］．設(shè)f1，f2是兩幅待拼接的圖像，將f1和f2在空間疊加，則融合后的圖像像素f 可表示為

設(shè)圖像間重疊區(qū)域?qū)挾葹長(zhǎng)，兩幅圖像重疊區(qū)x 軸和y 軸的最小和最大值分別為xmin、xmax和ymin、ymax，則過渡因子的計(jì)算方法為

式中:β，(1 -β)表示過渡因子，一般與重疊區(qū)域的寬度有關(guān)，且0 ＜β ＜1．在重疊區(qū)域中，β 由1 漸變到0，(1 -β)由0 漸變至1，由此實(shí)現(xiàn)了重疊區(qū)域中由f1到f2的平滑過渡．

4 實(shí)驗(yàn)結(jié)果及分析

圖7 ～9 是應(yīng)用本文算法和歸一化相似性算法(WC)分別對(duì)正常圖像、扭曲圖像及不同角度的旋轉(zhuǎn)圖像拼接效果進(jìn)行對(duì)比．從圖可看出，當(dāng)圖像沒有旋轉(zhuǎn)、扭曲或者色度均勻時(shí)不同算法的拼接效果相差很小;當(dāng)圖像有扭曲現(xiàn)象時(shí)本文算法對(duì)圖像的拼接處理更優(yōu);當(dāng)圖像旋轉(zhuǎn)角度達(dá)到30°時(shí)，NC 算法出現(xiàn)了圖像失配情況，而本文算法仍能夠?qū)崿F(xiàn)對(duì)圖像的拼接．實(shí)驗(yàn)結(jié)果表明:本文算法對(duì)圖像有更好的旋轉(zhuǎn)適應(yīng)性．

圖7 原圖像采用本文算法與NC 算法的拼接比較Fig．7 Original image using the algorithm of this paper compared with NC algorithm of joining together

圖8 扭曲30 度圖像采用本文算法與NC 算法的拼接比較Fig．8 Distorted 30° images using this algorithm and the algorithm of NC joining together

圖9 旋轉(zhuǎn)30 度圖像采用本文算法與NC 算法拼接比較Fig．9 30° rotation image compared with NC algorithm joining together using the algorithm of this paper

5 結(jié)論

提出了一種SPZMA 圖像拼接算法，該算法利用改進(jìn)的SURF 算子獲取圖像中的特征點(diǎn)，計(jì)算以特征點(diǎn)為中心鄰域窗口的偽Zernike 矩，獲得各個(gè)特征點(diǎn)鄰域偽Zernike 矩的Bray-Curtis 相似性度量的初始匹配點(diǎn)對(duì)，利用RANSAC 算法剔除偽特征點(diǎn)對(duì)，之后對(duì)輸入圖像作幾何變換進(jìn)行配準(zhǔn)，融合重疊區(qū)域，獲得良好的拼接圖像．實(shí)驗(yàn)結(jié)果表明:提取特征點(diǎn)均勻、準(zhǔn)確、迅速;圖像配準(zhǔn)時(shí)對(duì)平移、旋轉(zhuǎn)以及噪聲均具有魯棒性．

［1］ KNOPP J，PRASAD M，WILLEMS G，et al． Hough transform and 3D SURF for robust three dimensional classification［C］//Computer Vision-ECCV 2010-11th．Berlin Heidelberg:Springer，2010:589 -602．

［2］ WONG T C K． Box filter structure:U． S，Patent，7708，883［P］． 2010 -05 -04．

［3］ REININGHAUS J，KOTAVA N，GUNTHER D，et al．A scale space based persistence measure for critical points in 2d scalar fields［J］． IEEE Transactions on Visualization and Computer Graphics． 2011，17(12):2045 -2052

［4］ HEROUT A，HRADI? M，ZEM?íK P． EnMS:early non-maxima suppression［J］． Pattern Analysis and Applications，2012，15(2):121 -132．

［5］ 紀(jì)東升． 基于特征點(diǎn)的圖像拼接算法研究［D］． 蘭州:蘭州理工大學(xué)計(jì)算機(jī)與通信學(xué)院，2011．

［6］ SATTLER T，LEIBE B，KOBBELT L． Fast imagebased localization using direct 2D-to-3D matching［C］//IEEE International Conference on Computer Vision (ICCV)． Spain:IEEE Press，2011:667 -674．

［7］ YANG B，LI S． Multifocus image fusion and restoration with sparse representation［J］． IEEE Transactions on Instrumentation and Measurement，2010，59(4):884 -892．

［8］ ARIVAZHAGAR S，ANISHA J P． Image fusion using spatial unmixing［C］//Signal Processing Image Processing ＆ Pattern Recognition (ICSIPR)，2013 International Conference on IEEE． Coimbatore:IEEE Press，2013:238 -242．