李 壯,雷志輝,于起峰
國防科學技術大學航天與材料工程學院,湖南長沙410073
基于梯度徑向夾角直方圖的異源圖像匹配
李 壯,雷志輝,于起峰
國防科學技術大學航天與材料工程學院,湖南長沙410073
提出一種圖像全局特征描述方法——梯度徑向夾角金字塔直方圖。該特征不受圖像灰度級非線性變換和圖像旋轉的影響,具有較強的抗噪聲能力。基于該特征的匹配算法能夠對存在旋轉變換的異源圖像進行匹配。用SAR圖像和可見光圖像對匹配算法進行測試,結果表明基于梯度徑向夾角金字塔直方圖匹配方法的匹配成功率優于傳統方法。由于對基準圖進行離線預處理,同時避免角度空間搜索,算法的實時計算時間比傳統方法小3個數量級。
異源圖像;圖像匹配;旋轉不變;梯度矢量;徑向矢量;金字塔直方圖
由不同類型傳感器獲得的圖像被稱為異源圖像。由于不同傳感器間的成像特性存在差別,同一場景在異源圖像上可能呈現完全不同的圖像,因此傳統的同源圖像匹配方法一般無法直接應用于異源圖像。
國際上主流的異源圖像匹配方法可以分為兩大類,基于特征的匹配方法和基于區域的匹配方法。基于特征的匹配方法包括特征提取、特征描述、特征匹配三個步驟。理論上,基于特征的匹配方法對圖像變形有較好的適應性,但是在異源圖像中提取同名特征本身也是一個難題,目前的算法通常都只是針對某類圖像有效,而對于噪聲污染嚴重的圖像如SAR圖像等,尚無較好的同名特征提取方法。傳統的特征描述方法有小波特征[1]、形狀上下文[2]、矩方法[3]、SIFT[4]、PCASIFT[5]、SURF[6]等,文獻[7]中對不同的特征描述方法進行了分析和比較。這些特征描述方法只適用圖像描述區域較小,內容簡單,區域統計特性明顯的情況。例如SIFT特征的有效描述區域在幾十個像素內,且要求該區域內有較明顯的主方向。由于異源圖像提取的特征存在大量非同名點,因此傳統的特征匹配方法的有效性也難以保證。
基于區域的匹配方法通常將整幅圖像或圖像子區[8]的灰度或梯度信息進行匹配。此類方法中利用了圖像的所有的空間信息和灰度信息,抗噪性較基于特征的方法更高。常用的基于區域的匹配方法包括灰度互相關方法[9]、梯度相關方法[10]、互信息方法[11],這些方法的主要缺點是對圖像旋轉及灰度變化比較敏感,另外,由于計算相似度時需要用到全部像素信息,因此計算量也比較大,盡管有優化的匹配空間搜索策略[11-12],其優化后的時間仍難滿足實時性的要求。
基于區域的匹配方法能夠克服局部噪聲,而基于特征的匹配方法能夠適應圖像變形。綜合二者的優點,提出一種旋轉不變的圖像全局特征描述方法——梯度-徑向夾角金字塔直方圖(GRAPH),并將該特征應用在異源圖像匹配中。同時重點介紹梯度-徑向夾角金字塔直方圖的計算方法,在得到圖像的梯度-徑向夾角金字塔直方圖后,對基準圖不同窗口的特征向量與實時圖的特征向量進行比較,定義距離最近者為實時圖在基準圖中的最佳匹配。
建立梯度-徑向夾角金字塔直方圖分為下面五個步驟:選取有效區域;計算徑向矢量(radius direction vector,RDV);計算歸一化梯度矢量(normalized gradient vector,NGV);計算梯度-徑向夾角(gradient radius angle,GRA);用金字塔直方圖(pyramid histogram,PH)統計有效區域內梯度-徑向夾角的分布情況。
為使圖像特征描述具有旋轉不變性,首先需要選取旋轉不變區域作為計算特征描述的有效區域。在所有形狀中,只有圓形能夠滿足旋轉不變的性質,因此將特征描述的定義域設定為圓形。不失一般性,可以將圖像中心看作圖像旋轉的軸心,則取以圖像中心點為圓心的最大內接圓作為有效區域。如圖1,白色圓形區域為有效區域,圖像旋轉前后有效區域內部的像素能夠一一對應,而其他形狀如灰色矩形區域內的像素無法對應。定義有效區域內的像素為有效像素,下文中所有的計算都是針對有效像素進行,不再另加說明。

圖1 有效區域選取示意圖Fig.1 Select valid region
首先定義本文對矢量操作的標記符號。對于矢量v=a+bj,記矢量長度為矢量方向為∠v=acrtan(b/a)。對矢量旋轉操作記為Rotate(v,θ),其中θ為順時針旋轉角度。下文中vr表示徑向矢量,vd表示梯度矢量。
設圖像有效區域的中心點為O(x0,y0),則對于任意有效像素A=I(x,y),定義該點的徑向矢量為當圖像旋轉θ度后,A對應像素AR的徑向矢量為
文獻[10]中提出利用歸一化梯度矢量進行匹配的方法。圖像中一點 I(x,y)處的歸一化梯度矢量定義為

梯度矢量歸一化操作本質上是去除梯度矢量的強度信息,只保留梯度矢量方向。文獻[10]主要用歸一化梯度矢量解決同源圖像存在亮度變化時的匹配問題,并沒有考慮異源圖像灰度反轉或存在密集噪聲的問題,另外圖像中存在梯度方向不穩定像素,如湖泊,其內部梯度很小但是通常不為零,則梯度方向一般僅受噪聲影響,分布非常雜亂。將這類不穩定像素的歸一化梯度矢量加入到圖像匹配算法中,必將影響匹配的可靠性。
為了解決異源圖像中灰度反轉,即同名點的梯度方向可能相反的問題,一種直觀的想法是將圖像梯度由[0,2π)規定到[0,π),本文中將方向規定化操作放在計算梯度-徑向夾角后進行,具體見2.4節公式(7)。
為了消除噪聲影響,可以采用低通濾波。實際上,濾波操作還會增大梯度方向一致區域,利于圖像匹配。
為了消除不穩定像素的影響,設定閾值 Th,舍棄梯度強度小于閾值的像素點。為了使閾值適用于不同類型圖像,先將圖像進行濾波和灰度線性增強處理。在灰度線性增強中,去除圖像灰度直方圖中最亮和最暗的像素各1%,將剩余的灰度級線性拉伸到[0,255]區間。梯度強度閾值 Th設為20,對于不同類型圖像都得到較好結果。
圖像旋轉θ度后,I(x,y)對應的像素點為IR(x1,y1),其歸一化梯度矢量也將旋轉θ度。

定義vr與vd的夾角α為梯度-徑向夾角。圖像旋轉θ度后,與A對應的像素AR的梯度-徑向夾角為

由上式可知,梯度-徑向夾角不受圖像旋轉影響。已知歸一化梯度矢量和徑向矢量,則可通過反余弦函數計算梯度-徑向夾角。

規定梯度矢量相位超前徑向矢量相位時,梯度-徑向夾角為正,將α量化到(-π,π),記為α2π,有如下計算公式,其中vdx、vdy分別為vd的實部和虛部,vry、vrx為 vr的實部和虛部。vdx×vry-vdy× vrx>0表示vd相位超前vr,vdx×vry-vdy×vrx<0表示vd相位滯后vr。

為了適應異源圖像中梯度反向的情況,需將梯度-徑向夾角量化到[0,π)。

需要特別說明,任意一點的徑向矢量只與該點在窗口中的幾何位置有關,故可以預先將窗口內每一點的歸一化徑向矢量制成二維查找表。在對基準圖進行匹配搜索時,每個窗口只需查表就可確定其中每點的歸一化徑向矢量。圖像匹配時,同樣可以預先將反余弦函數制成查找表,通過查表加快算法速度。
圖像旋轉后,由于對應點的空間坐標發生改變,無法直接比較對應像素的梯度-徑向夾角。一種合理的方法是采用直方圖表示圖像的區域統計特征。建立梯度-徑向夾角的直方圖,則該直方圖具有旋轉不變特性。由公式(6)、(7)將夾角量化到[0,π)區間。本文中將角度空間劃分12個區間,每個區間對應直方圖一個柵格,則可以用12維向量表示直方圖。
實際應用時,僅用12維向量表示一幅圖像,則向量的分辨能力有限,無法實現可靠匹配。文獻[13]在研究圖像分類時,提出一種梯度方向的金字塔直方圖特征,能夠通過金字塔層數的改變調整對圖像的表示能力。為了增加向量維數,提高特征的表示能力,同時保證特征描述的旋轉不變性,本文采用一種圓形金字塔分層方法。圖2是3層圓環金字塔示意圖,其中L1層每個區域面積為L0層區域面積的一半,L2層每個區域面積為L1層區域面積的一半,依此類推。分別在L0、L1、L2層每個子區計算梯度-徑向夾角直方圖,由于不同層上的區域直方圖統計的區域面積(有效像素)不同,用直方圖向量除以該統計區域面積得到歸一化直方圖(為了計算時可以用整數表示直方圖向量,再乘上L0層區域面積)。將全部歸一化直方圖向量首尾拼接,得到圓環區域金字塔直方圖。圖3為5層圓形金字塔梯度-徑向夾角直方圖特征,由372維向量表示。

圖2 圓形金字塔直方圖Fig.2 Round pyramid histogram

圖3 梯度-徑向夾角金字塔直方圖特征Fig.3 GRAPH descriptor
理論上,梯度-徑向夾角金字塔直方圖特征是一種真正的圖像旋轉不變稀疏特征描述,它不依賴于主方向的選取,比傳統的方法有更強的適應性。梯度-徑向夾角金字塔直方圖特征能夠實現圖像灰度非線性變化,梯度方向反向等復雜條件下同源或異源圖像間的匹配。
在圖像匹配的許多應用中,只需要得到一幅圖像中心點在另一幅圖像中的對應位置,無需獲得兩幅圖像的全部變換參數下精確配準。這種情況下,可以使用中心點匹配模型。中心點匹配模型中,有實時圖 I1和基準圖 I2,通常實時圖場景完全包含在基準圖場景中,匹配目的是要找出實時圖中心在基準圖中的對應位置。以圖像中心為原點,用 TRT={x,y,θ}表示對基準圖的平移-旋轉變換(RT變換),IRT表示全體RT變換的集合,則圖像匹配的目的就是尋找使某種條件達到極值的 TRTm。即

S upport(I1,TRT(I2))為對基準圖做 RT變換 TRT時,對匹配的一種支持度。
如果能夠找到一種旋轉不變的圖像描述D,使得S upport(D(I1),D(TRT(I2)))=S upport (D(I1),D(TT(I2))),其中 TT={x,y}表示對基準圖的平移變換,則只需進行平移空間的搜索匹配,就能夠得到最佳匹配位置。梯度-徑向夾角金字塔直方圖正是這樣一種旋轉不變圖像描述。對于圖像 I,其梯度-徑向夾角金字塔直方圖特征D(I)是一維向量。則可以設定 S upport(D(I1), D(TT(I2)))=-Dis(D(I1),D(TT(I2))),其中Dis(·)是向量間的χ2距離[14]。對于向量 va= {ai},vb={bi},有

圖4中給出基于梯度-徑向夾角金字塔直方圖的圖像匹配算法,該算法適用于同源圖像或者異源圖像間的匹配。低通濾波和灰度線性增強是非常重要的兩個預處理過程,試驗中發現,對于SAR圖像和可見光圖像的匹配,如果不做這兩步,將會影響匹配的效果。雙查找表的加入大大提高了算法運算速度。整個匹配過程采用平移位置上搜索匹配的策略,搜索步長設為2像素。算法中的灰色區域表示該步驟可以在任務開始前離線處理和加載。因此,實時處理的步驟只有實時圖的處理和直方圖距離計算,大大減少了實時任務的計算壓力。

圖4 圖像匹配算法流程圖Fig.4 Image matching algorithm flow chart
理論上,梯度-徑向夾角金字塔直方圖不受圖像旋轉和灰度級非線性變換的影響,為了驗證這兩個性質,下面用仿真圖測試梯度-徑向夾角金字塔直方圖特征的抗旋轉和灰度變化性能。圖5中,第一行是輸入圖像和對應的梯度-徑向夾角金字塔直方圖特征。第二行是對輸入圖像旋轉并作非線性灰度級變換和對應的梯度-徑向夾角金字塔直方圖特征。圖中兩個直方圖相似度很高,通過計算直方圖距離,發現它們之間存在細微差別。存在差別的主要原因是計算梯度時采用了3×3方形窗口,在圖像旋轉時,對應像素點的梯度矢量除旋轉外還會存在由于方形窗口引入的誤差,考慮到此差異很小,幾乎不會對匹配造成影響,因此并沒有采用更復雜的方法,如大尺度圓形窗口來消除它。圖5中的試驗說明梯度-徑向夾角金字塔直方圖特征能夠克服圖像灰度級非線性變換和圖像旋轉的影響。

圖5 圖像旋轉和灰度級變換下的梯度-徑向夾角金字塔直方圖特征Fig.5 GRAPH feature after image rotation and grey level transform
為了測試梯度-徑向夾角金字塔直方圖特征對不同圖像的表達力和分辨能力,對多幅衛星圖像計算特征進行描述。采用χ2距離計算不同圖像的梯度-徑向夾角金字塔直方圖特征間的差異程度,見表1,表內數值為向量距離,無量綱。試驗說明,對于任意兩幅不相關的圖像,它們的直方圖距離比較大,對于任意選擇的7幅圖像,直方圖距離最小的為6 782。
對一幅衛星圖像及其旋轉30°和60°后的圖像計算梯度-徑向夾角金字塔直方圖特征描述如圖6。采用χ2距離計算得到直方圖間距離如表2,其中最大的距離為648。與表1中的距離相比較,圖像旋轉產生的向量距離小一個量級以上。

表1 不同場景圖像的直方圖特征距離Tab.1 Distances of different GRAPHfeature

圖6 圖像旋轉不同角度的梯度-徑向夾角金字塔直方圖特征描述Fig.6 GRAPH feature descriptor respected to different rotate angles

表2 圖像旋轉不同角度的梯度-徑向夾角金字塔直方圖特征距離Tab.2 Distances of different GRAPHfeatures respected to different rotate angles
基于梯度-徑向夾角金字塔直方圖特征進行圖像匹配,首先測試匹配算法對同源圖像的效果。基準圖為一幅遙感照片,實時圖是從基準圖中選取的一小塊區域,并進行了角度旋轉。匹配結果如圖7(a),實時圖和基準圖間只有平移和旋轉變化,白色圓環的中心對應實時圖中心在基準圖中的匹配位置。圖7(b)是過匹配點沿 x軸方向各窗口與實時圖直方圖相似性(距離倒數)曲線。匹配真值為(406,241),梯度-徑向夾角直方圖特征匹配結果為(406,240)。圖7的試驗驗證了本文方法原理上的正確性,即對同源圖像能夠得到正確結果。


圖7 旋轉圖像匹配Fig.7 Matching rotated images
為了測試梯度-徑向夾角金字塔直方圖特征對異源圖像的匹配效果,選取美國Sandia國家實驗室提供的SAR圖像做實時圖,Google Earth中的光學衛星圖像做基準圖,兩圖存在平移旋轉變化,另外還有一定程度的圖像畸變。為了考核匹配結果是否正確,手工選擇多組同名特征點,用最小二乘方法得到最優匹配解作為匹配點的參考值。由于手選同名點的存在誤差,故無法說明這兩個結果哪個更真,但是只要結果接近,一般坐標差異在5個像素內,即可認為匹配結果可以接受,本方法只考核匹配正確性,并不考核匹配精度。采用梯度-徑向夾角金字塔直方圖特征匹配結果如圖8(a)、9(a)。中心點匹配模型,圖8(a)中匹配參考值為(363,220),本文方法匹配結果為(361,221),匹配結果與參考值距離約2個像素,故匹配結果正確。過匹配點沿 x軸方向各窗口與實時圖直方圖相似性曲線如圖8(b),圖8(c)、(d)分別為互信息和文獻[8]方法匹配結果,匹配位置錯誤。圖9(a)中匹配參考值為(116,330),本文方法匹配結果為(116,331),匹配結果正確。過匹配點沿 x軸方向各窗口與實時圖直方圖相似性曲線如圖9(b),圖9(c)、(d)分別為互信息和文獻[8]方法匹配結果,匹配位置錯誤。


圖8 異源圖像匹配Fig.8 Matching multi-sensor images
圖8、圖9說明本文方法能夠適用于異源圖像的匹配,但是相似曲線主次峰比值較小,表明異源圖像匹配結果可靠性沒有同源圖像匹配可靠性高。準備了一組待匹配圖像,包括30對SAR圖像和光學圖像,SAR圖像作為實時圖,尺寸為300× 300像素,光學圖像作為基準圖,尺寸為600×600像素。分別對其中的SAR圖像旋轉15°、30°、45°得到另外三組待匹配圖像。將本文方法與互信息及文獻[8]方法進行了對比。本文方法將基準圖直方圖計算作為離線處理,不計入計算時間。互信息方法和文獻[8]方法分別進行平移搜索和平移加旋轉搜索,其中旋轉搜索步長為10°。匹配成功率為成功匹配圖像數比總的圖像數。三種方法的匹配成功率和計算時間見表3。在圖像只存在平移變換時,文獻[8]方法只進行平移搜索時的成功率最高,本文方法次之,而文獻[8]方法進行平移加旋轉搜索時,其成功率大幅降低,低于本文方法。對旋轉15°、30°、45°的圖像進行匹配時,本文方法成功率遠高于互信息方法和文獻[8]方法。實時計算時間上本文方法遠低于了互信息方法和文獻[8]方法,降低了2個數量級以上。這是因為本文方法無需旋轉搜索,且將基準圖進行了離線計算,實時計算時只需獲得實時圖的梯度徑向夾角直方圖,并進行直方圖比較操作。而另外兩種方法對每個窗口計算時都同時需要實時圖和對應基準圖窗口內的全部像素信息,無法進行離線預處理。

圖9 異源圖像匹配Fig.9 Matching multi-sensor images

表3 不同匹配方法對比Tab.3 Comparison of different matching methods
上面試驗中的異源圖像之間存在復雜的幾何形變,憑經驗給出的匹配真值只能用以判斷匹配正確與否,不能用來考核算法匹配精度。為此,制作仿真異源圖對算法精度進行考核。首先選擇內容豐富的光學衛星圖像作為基準圖。取基準圖上的一部分區域,進行非線性灰度變換,加入斑點噪聲并隨機旋轉一定角度,得到仿真實時圖。通過對多組仿真圖像的匹配試驗得到匹配誤差均值為0.9像素,誤差方差為0.8。
本文提出了一種基于梯度-徑向夾角金字塔直方圖特征的旋轉不變圖像描述方法,并將它應用于任意旋轉角度下的異源圖像匹配。試驗說明了該特征描述方法能夠克服圖像旋轉和灰度變化。在存在圖像旋轉、平移和一定程度圖像畸變的異源圖像上進行了測試,結果表明該方法能夠實現大部分異源圖像間的正確匹配。另外本文方法在處理時間上大大優于各種常有的匹配方法。只利用梯度-徑向夾角金字塔直方圖特征進行匹配,尚無法解決圖像存在較大尺度變化和形狀畸變的情況下的匹配。下一步研究工作重點是提高匹配成功率,同時爭取解決異源圖像同時存在尺度和旋轉變化時的匹配問題。
致謝:感謝美國Sandia國家實驗室提供SAR圖像,感謝Google提供光學衛星圖像。
[1] CHEN Ying,YE Qin,ZHONG Zhiyong.Research on the Matching between Radar Image and Optical One Based on Wavelet[J].Acta Geodaetica etCartographica Sinica, 2000,29(3):245-249.(陳鷹,葉勤,鐘智勇.基于小波變換的雷達與光學影像匹配算法研究[J].測繪學報,2000, 29(3):245-249.)
[2] BELONGIE S,MALIK J,PUZICHA J.Shape Matching and Object Recognition Using Shape Contexts[J].IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002,24(4):509-522.
[3] VAN GOOL L,MOONS T,UNGUREANU D.Affine/ Photometric Invariants for Planar Intensity Patterns[C]∥Proceedings of the 4th European Conference on Computer Vision(ECCV′96).Cambridge:Springer,1996:642-651.
[4] LOWE D G.Distinctive Image Features from Scale-invariant Keypoints[J].International Journal of Computer Vision, 2004,60(2):91-110.
[5] KE Y,SU KTHANKAR R.PCA-SIFT:a More Distinctive Representation for Local Image Descriptors[C]∥Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition(CVPR′04).Washington: IEEE Computer Society,2004:506-513.
[6] BA Y H,TU YTELAARS T,VAN GOOL L.SURF: Speeded Up Robust Features[C]∥Proceedings of the 9th European Conference on Computer Vision(ECCV′06). Graz:Springer,2006:404-417.
[7] MIKOLAJCZYK K,SCHMID C.A Performance Evaluation of Local Descriptors[C]∥Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition(CVPR′03).Madison:IEEE Computer Society, 2003:257-263.
[8] ZHANG Shaoming,CHEN Yingying,LIN Yi.A Robust Matching Algorithm for SAR Image with Multiple Subareas[J].Acta Geodaetica et Cartographica Sinica,2007, 36(4):406-413.(張紹明,陳映鷹,林怡.用于末制導的SAR圖像多子區實時匹配算法[J].測繪學報,2007,36 (4):406-413.)
[9] YU Qifeng,SHANG Yang.Videometrics:Principles and Researches[M].Beijing:Science Press,2009:96-102.(于起峰,尚陽.攝影測量學原理與應用研究 [M].北京:科學出版社,2009:96-102.)
[10] FITCH A J,KADYROV A,CHRISTMAS W J,et al. Orientation Correlation[EB/OL].[2010-01-10].http:∥citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1. 64.5662&rep=rep1&type=pdf.
[11] MAES F,COLLIGNON A,VANDERMEULEN D,et al.Multimodality Image Registration by Maximization of Mutual Information[J].IEEE Transactions on Medical Imaging,1997,16(2):187-198.
[12] LAMPERT C H,BLASCH KO M B,HOFMANN T.Beyond SlidingWindows:ObjectLocalizationby Efficient Subwindow Search[C]∥Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition(CVPR′08).Los Alamitos:IEEE Computer Society,2008:1-8.
[13] BOSCH A,ZISSERMAN A,MUNOZ X.Representing Shape with a Spatial Pyramid Kernel[C]∥Proceedings of the 6th ACM International Conference on Image and Video Retrieval(CIVR′07).Amsterdam:[s.n.],2007: 401-408.
[14] MARTIN D,FOWL KES C,MAL IK J.Learning to Detect Natural Image Boundaries Using Local Brightness, Color and Texture Cues[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence,2004,26(5): 530-549.
Matching Multi-sensor Images Based on Gradient Radius Angle Pyramid Histogram
LI Zhuang,LEI Zhihui,YU Qifeng
College of Aerospace and Material Engineering,National University of Defense Technology,Changsha 410073,China
A novel image feature named gradient radius angle pyramid histogram(GRAPH)is proposed.The GRAPH feature is invariable to nonlinear transform of grey level and image rotation,and is robust to noise.Based on GRAPH feature,we get a robust and fast image matching algorithm that can match rotated and shifted multi-sensor images.SAR images and optical images are used to test this matching method.The results show that the new matching algorithm has higher success rate than traditional method.As the reference image can be processed offline,the processing time is less then one in a thousand of the traditional used time.
multi-sensor image;image matching;rotate-invariant;gradient vector;radial direction vector; pyramid histogram
LI Zhuang(1982—),male,PhD candidate,majors in image registration and target recognition.
1001-1595(2011)03-0318-08
TP391
A
國家863計劃(2007AA12Z121)
(責任編輯:雷秀麗)
2010-02-24
2010-07-31
李壯(1982—),男,博士生,主要研究方向為圖像匹配和目標識別。
E-mail:lizhuang2007@hotmail.com