吳偉平,金龍旭,閆得杰,王 棟
(1.中國科學院 長春光學精密機械與物理研究所,吉林 長春130033,2.中國科學院大學,北京100039)
圖像配準技術主要可以分為基于灰度的方法和基于特征的方法兩大類。后者由于提取特征后僅對特征進行計算,相對前者其計算量較少,對噪聲、光照、視角和尺度變化不敏感,算法效率及配準精度高,具有很好的魯棒性,因而成為當前圖像配準領域的主要研究方向。早期基于特征的算法包括Morvec、Harris[1]等方法。
Bay等 提 出 了 著 名 的SURF (speeded up robust features)[2]算法,SURF受滑動角度的影響,當配準圖像間的旋轉方向是滑動角度的整數倍時,SURF的表現略有下降。為改進這一問題,本文通過優化不變矩系數矩陣提高了描述子方向計算精度,分別設計了兩種基于圖像亮度質心不變矩的描述子方向計算方法,并將其與SURF 描述子方向的計算精度和效率做了詳細地比較。
為保證計算精度,描述子方向的計算一般采用以特征點為 中 心 的 圓 形 圖 像 區 域。E.Rublee 等 在 設 計ORB[4](oFAST and rBRIEF)描述子中采用了Rosin對圖像矩的定義:mpq=∑xpyqI(x,y),其中p、q取0或1,其亮度質心為:c=(m10/m00,m01/m00),從特征點o到亮度質心c建立一個代表圖像片方向的向量oc,則向量oc與X 軸的夾角可以表示為:c=arctan (m01/m10)。當m10接近0時,該方法將不具備穩定性。因此,當m10小于文中設定的閥值時,該特征點不具備方向性,將該點丟棄不參與配準,在實驗過程中這種情況是非常罕見的。由θ的定義可知θ∈ (-π/2,π/2),而特征點描述子的方向區間需要達到 (-π,π)。由于特征點是極值點,因此m00=∑xyI(x,y)>0。可以根據m01的符號將θ擴展到 (-π,π),即可分為如下3種情況計算

在ORB配準方法中,圖像片直接選取特征點(xc,yc)為圓心,半徑為r以內的所有完整的像素點 ((x,y)|(x-xc)2+ (y-yc)2≤r2且x,y ∈R)進行圖像矩計算,而忽略了像素大小和圓周邊緣上的點對圖像矩的影響。這種忽略方式簡化了計算但損失了較多方向精度。文中對這種方向矢量計算方法進行了改進,設計了軸向密集插值法和面積積分法。
當圖像旋轉不同角度進行亮度質心矩計算時,圖像矩的定義可知圖像圓周上點的權重比內部像素點的權重大,對圖像矩的影響也最大,因此圓周上像素是圖像片中最需要精確計算的部分。圖像片圓周上的點可通過插值的方法得到[5],這些點的像素值不僅和圓周內部的點有關,還和圓周外相鄰近的點有關。軸向密集插值法正是基于這種突出圓周像素的思想而設計的。將整個圓周從0°開始每間隔π/4分為一個區域 (如圖1所示),(π/4,π/2)圓弧所對應X 軸的長度為X0,對應Y 軸的長度為Y0,由于X0>Y0,在相同插值間隔條件下,選擇插值點更密集的軸向可以使插值結果更均勻。因此采用X 軸方向插值計算,先用式(1)求出所有插值點橫坐標,再用式 (2)求得相應點的縱坐標;同理,在 (0,π/4)圓弧上則采用Y 軸方向插值計算,按照這種方式將整個圓周進行分段插值,得到圓周上的全部插值點


圖1 半徑為9的軸向密集插值法
則圓弧上的插值點亮度用軸向線性插值的方法計算,如式 (3)所示


由于軸向密集插值法的計算點在圓域上分布并不絕對均勻,且該方法將像素點看做有亮度無大小的點。這樣做在計算時很容易區分該點是否在圓域內外,但其本質并不符合成像原理。數字圖像中的像素代表了一定面積物體的總亮度,是有大小的。因此本文設計了考慮像素大小的面積積分法來計算描述子方向。面積積分法將每個像素看成1*1的正方形,并假定每個像素的亮度在其正方形域內均勻分布。本方法將方向矢量計算的圓域中心設為坐標軸原點,則特征點坐標為x=-0.5,x=0.5,y=0.5,y=-0.5四條直線圍成的面積。如圖2所示,設點 (x,y)的積分系數為Qxy,代表該點在圓域中所占的面積,則面積積分法圖像矩公式變為mpq=∑QxyxpyqI(x,y)。

圖2 面積積分法陰影區
每個像素對應的系數為該像素在圓域中的面積。其中圓域內非陰影區方格的系數為1,其坐標 (x,y)與該方格代表的像素坐標一致。而陰影區方格的系數為該陰影區的面積,坐標為該陰影區的亮度質心的坐標。陰影區的面積和坐標均在配準前計算得到,然后將相應項的乘積組成系數矩陣保存在配準程序中,因此不占用配準時間。如表1所示半徑為4的圓域所對應的面積積分法系數矩陣,該矩陣采用定積分的方法計算得到。

表1 半徑為4的面積積分法系數矩陣
陰影區坐標 (xi,yi)采用幾何形心的計算方法,先利用分割法將陰影區分割為矩形和弓形,采用定積分得到弓形面積,再加上矩形面積得到陰影區面積。以圖3中陰影區域的xi坐標的求解過程為例,應用式 (4),其中dS表示積分變量dx 所對應的陰影區面積,S 為整個陰影區的面積,積分精度選擇10-7。

圖3 陰影區坐標求解

配準的硬件結構如圖4所示,實驗采用2臺DALSA 公司的Falcon 4M60彩色面陣相機拍攝圖像,采用加拿大IO Industries公司的CLSAS圖像處理卡同時采集兩臺相機輸出的圖像。該采集卡可同時采集兩個數據流,最大采集速度可達到450 MBps,并將采集的數據直接存儲到4TB 的SAS高速磁盤陣列中,可通過程序選擇部分圖像數據通過PCI總線上傳至計算機進行配準處理。配準使用的計算機采用Intel 4960X,3.6GHz的CPU、4GB 內存和XP SP3的操作系統,配準軟件采用C++語言編寫[6-8]。

圖4 配準硬件結構
在配準算法的實現中,首先構建積分圖和金字塔,采用Hessian矩陣濾波器對圖像進行特征點檢測,Hessian矩陣H(x,σ)定義為

式中:σ為當前尺度,Lxx(x,σ)是高斯二階偏導數與圖像I在點X 處的卷積,Lxy(x,σ)和Lyy(x,σ)也以此類推。構建5組不同尺度的金字塔數據,每組由4層構成,各層由大小不同箱式濾波器近似Hessian矩陣濾波器,具體濾波器大小如表2所示。這樣保證了除第一組的第一層和第五組的第三、四層以外,每層圖像剛好既作為中間層有作為邊緣層參與特征點檢測計算,這樣可以減少重復計算次數,提高配準效率。然后對每組的第二、三層進行3*3*3鄰域的非極大值抑制得到候選點。利用閥值法根據候選點二階梯度篩選得到特征點集合。采用半徑為9σ的軸向密集法和面積積分法計算特征點的方向。在不同尺度σ下計算計算描述子方向時,將以σ為邊長的正方形作為1個偽像素,正方形內全部像素的平均值作為該偽像素的亮度值進行方向計算,這種方式恰好可以利用SURF 的積分圖加快計算速度。

表2 金字塔構建參數
得到描述子方向后,以特征點為中心,描述子方向為主方向,選取邊長為20σ 的正方形區域,并將其劃分為4*4個子區域,對每個子區域計算x、y方向的Haar小波響應dx和dy(小波模板大小為2σ*2σ),并累加得到4維向量V(∑dx,∑dx,∑|dx|,∑|dx|),所有的子區域連接起來得到64維的特征描述子進行配準。
在實時配準實驗中,對兩臺視角基本相同的相機采集的圖像進行配準實驗,相機2在0-90°區間每隔15°旋轉一次[9,10]。整個實驗過程中,采用軸向密集插值法和面積積分法計算描述子方向矢量的配準過程可穩定高效地進行配準工作。配準 過 程 采 用RANSAC (random sample consensus)方法隨機抽樣一致方法對配準結果進行篩選,并根據RANSAC方法計算出的內點數計算配準的可重復率。
面積積分法配準實驗效果如圖5所示。

圖5 面積積分法配準實驗效果
設在基準圖像和待配準圖像中探測到的特征點個數分別為n1和n2。n=min (n1,n2)。m 是被證實匹配成功的點對數,即文中采用RANSAC 方法計算出的內點數[11]。則可重復率R=m/n。計算n個方向矢量的時間為t。表3中是對相機1與相機2所采集圖像的配準結果。實驗表明,采用面積積分法和軸向密集插值法的方向矢量計算時間大約是SURF方法的1/4,其配準可重復率方面,如圖6 所示,軸向密集插值法在45°角上的表現為3 種方法中最好的,而其它一些角度的表現則為三者中最差的,反映出了該方法在整個圓周上的表現不夠均衡,參數仍有很大的優化空間。面積積分法在各個角度上的表現具有很強的一致性,且平均重復率與SURF 法相當,較好地克服了SURF法在30°和60°的性能下降,并且獲得了較高的計算效率。

表3 實驗各旋轉角度配準結果

圖6 不同方向的配準重復率比較
為了檢測文中設計的方向計算方法在抗噪聲、視角變化、圖像壓縮等方面的表現,通過測試標準庫中的Trees、Leuven、和Ubc等圖像,表明采用亮度質心不變矩計算描述子方向在上述干擾條件下的總體表現可以達到SURF 算法的水平。其中,面積積分法的精度和穩定度都要好于軸向密集插值法,對視角變化和圖像壓縮的表現略好于SURF,而對強噪聲干擾方面略遜于SURF。
文中比較了兩種基于亮度質心不變矩的描述子方向改進計算方法和SURF方法的配準性能,兩種不變矩配準方法的計算速度均達到SURF 方法的4倍。在角度計算精度方面,面積積分法在整個圓周的表現要比SURF 和軸向密集插值法更加均衡。文中采用的系數矩陣的計算方式,在具有浮點型計算精度的同時,極大地優化計算速度,使得很多復雜的計算可以再配準前完成,不占用配準時間。文中軸向密集插值法的不均衡性說明系數矩陣仍有提升的空間,是下一步需優化的主要問題。
[1]LV Xuan,DUAN Huichuan.Multimodality medical image registration by mutual information and Harris corner detector[J].Computer Engineering and Design,2008,29 (4):998-1000 (in Chinese).[呂煊,段會川.基于Harris角點和最大互信息的多模醫學圖像配準 [J].計算機工程與設計,2008,29 (4):998-1000.]
[2]Bay H,Ess A,Tuytelaars T,et al.Gool,SURF:Speeded up robust features [J].Computer Vision and Image Understanding,2008,10 (3):346-359.
[3]Morel JM,Yu G.Asift:A new framework for fully affine invariant image comparison [J].SIAM Journal on Imaging Sciences,2009,2 (2)438-469.
[4]Rublee E,Rabaud V,Konolige K,et al.ORB:An efficient alternative to SIFT or SURF [C]//IEEE International Conference on Computer Vision,2011:2564-2571.
[5]FU Xiang,GUO Baolong.Overview of image interpolation technology[J].Computer Engineering and Design,2009,30(1):141-144 (in Chinese).[符祥,郭寶龍.圖像插值技術綜述 [J].計算機工程與設計,2009,30 (1):141-144.]
[6]MIAO Ligang.Image mosaicking and compositing algorithm for video surveillance [J].Chinese Journal of Scientific Instrument,2009 (4):857-861 (in Chinese).[苗立剛.視頻監控中的圖像拼接與合成算法研究 [J].儀器儀表學報,2009(4):857-861.]
[7]PENG Bo,HE Bin.Application and realization of FPGA in video mosaicking [J].Computer Engineering and Design,2013,34 (5):1635-1639 (in Chinese). [彭勃,何賓.FPGA 在視頻拼接中的應用與實現 [J].計算機工程與設計,2013,34 (5):1635-1639.]
[8]LI Jian,KONG Lingyin.Research and implementation of stereo match algorithm based on OpenCV [J].Computer Engineering and Design,2013,34 (2):566-569 (in Chinese).[李健,孔令寅.基于OpenCV 的立體匹配算法的研究與實現[J].計算機工程與設計,2013,34 (2):566-569.]
[9]LIU Tao,YU Zhongqing,MA Qianli.Target distance calculation based on parallel binocular vision [J].Qingdao University (Natural Science),2009,22 (1):59-62 (in Chinese).[劉濤,于忠清,馬千里.基于平行雙目視覺的目標距離計算[J].青島大學學報 (自然科學版),2009,22 (1):59-62.]
[10]TU J,TAO H,HUANG T.Online updating appearance generative mixture model for mean-shift tracking [J].Machine Vision and Applications,2009,20 (3):163-173.
[11]ZHOU Jianjun,OUYANG Ning,ZHANG Tong,et al.Image mosaic method based on RANSAC [J].Computer Engineering and Design,2009,30 (24):5692-5694 (in Chinese).[周劍軍,歐陽寧,張彤,等.基于RANSAC 的圖像拼接方法 [J].計算機工程與設計,2009,30 (24):5692-5694.]