余春超,楊智雄,夏宗澤,袁小春,嚴(yán) 敏
?
采用GPU并行架構(gòu)的基于互信息和粒子群算法的異源圖像配準(zhǔn)
余春超1,楊智雄1,夏宗澤2,袁小春1,嚴(yán) 敏1
(1.昆明物理研究所,云南 昆明 650223;2.國(guó)網(wǎng)遼陽(yáng)供電公司,遼寧 遼陽(yáng) 111000)
對(duì)于異源可見(jiàn)光與紅外圖像配準(zhǔn),以互信息為相似性度量條件,以粒子群算法為搜索算法,在搜索空間中搜索極大值點(diǎn),通過(guò)改變圖像分辨率,由粗到精,逐步實(shí)現(xiàn)可見(jiàn)光圖像與紅外圖像的配準(zhǔn),對(duì)粒子群算法、統(tǒng)計(jì)互信息和仿射變換3個(gè)部分的內(nèi)在并行性,利用CUDA C語(yǔ)言實(shí)現(xiàn)GPU-CPU異構(gòu)編程。實(shí)驗(yàn)證明,在不降低精度的前提下提高了算法效率,得到了很好的加速比,算法正確匹配率高,魯棒性好,計(jì)算效率高。
異源圖像;互信息;粒子群;圖像配準(zhǔn);GPU
異源圖像配準(zhǔn)[1]主要指不同類(lèi)傳感器圖像之間、或不同波段圖像之間的配準(zhǔn),如紅外與可見(jiàn)光圖像[2]、SAR圖像與光學(xué)圖像,以及醫(yī)學(xué)圖像中的CT與PET等。由于異源圖像成像機(jī)理不同,同一目標(biāo)在相同的外在條件下也會(huì)產(chǎn)生不同的特征表現(xiàn)形式,如圖1所示。
圖1所示為同一場(chǎng)景可見(jiàn)光(左)與紅外(右)圖像對(duì)比圖,由于紅外線的穿透性,在有較多煙霧遮擋的情況下紅外圖像可以捕捉到煙霧后面的影像信息;對(duì)于沒(méi)有被煙霧遮擋的物體,如圖像右側(cè)的房屋和前景地面,在可見(jiàn)光圖像(左)中表現(xiàn)得更為清晰,即邊緣梯度較大,而在紅外圖像(右)中則顯得較為模糊;此外,一些輻射能量較高的物體在紅外圖像中可能更為明顯[3]。
由于成像機(jī)理不同,局部相同信息過(guò)少,不能采用基于局部特征的配準(zhǔn)算法對(duì)異源圖像進(jìn)行配準(zhǔn)[4]。例如,基于特征點(diǎn)的算法,很難找到同名點(diǎn),或即便找到,其特征支撐區(qū)域也很可能不同;對(duì)于基于邊緣的算法,不同源的圖像邊緣很可能在位置或長(zhǎng)度上不一致,從而導(dǎo)致配準(zhǔn)失敗;基于局部區(qū)域的算法同樣包含上述問(wèn)題。
圖1 同一場(chǎng)景可見(jiàn)光(左)與紅外(右)圖像對(duì)比圖
針對(duì)異源圖像,本文采用基于全局灰度的配準(zhǔn)算法,從整體上搜索兩幅異源圖像的相似特征,避免了基于局部特征算法的不足。這類(lèi)算法的基本框架包含4個(gè)方面,分別為:特征空間、搜索空間、相似性度量以及搜索策略。
本文選取全局灰度信息作為算法的特征空間,即考慮圖像中所有像素的灰度值,且無(wú)需對(duì)圖像進(jìn)行預(yù)處理。由于考慮了所有灰度值,缺點(diǎn)是計(jì)算量較大,需要采用相關(guān)的優(yōu)化技術(shù)進(jìn)行優(yōu)化,本文在算法上采用粒子群搜索算法加速搜素,在硬件上利用GPU[5]的并行處理[6]能力對(duì)算法進(jìn)行并行設(shè)計(jì)。
針對(duì)本文算法,這里的搜索空間指配準(zhǔn)圖像在一定范圍內(nèi)進(jìn)行平移、旋轉(zhuǎn)、縮放后形成的一組圖像。建立搜索空間的目的是在其中利用搜索算法找出與待配準(zhǔn)圖像具有最大相似度的圖像,那么此圖像的變換參數(shù)即要求取的用于配準(zhǔn)的參數(shù),一般包括方向平移像素?cái)?shù)、方向平移像素?cái)?shù)、旋轉(zhuǎn)角度、縮放倍數(shù)。
互信息(Mutual Information)[7]是信息論里的一種度量,反映了被測(cè)事件集合之間的相關(guān)程度的大小[8]?;バ畔⒌闹翟礁撸硎臼录g的相關(guān)性也就越大。一般而言,信道中總是存在著噪聲干擾,信源發(fā)出的信號(hào)在信道中受到一定程度和形式的干擾,到達(dá)信宿后,由原來(lái)的變?yōu)?。通過(guò)后驗(yàn)概率(|),可以推測(cè)出原始信號(hào)的概率,而發(fā)出的先驗(yàn)概率為()。與的互信息可被定義為的后驗(yàn)概率與先驗(yàn)概率的比值取對(duì)數(shù)。在異源圖像配準(zhǔn)中,配準(zhǔn)圖像可以當(dāng)作,待配準(zhǔn)圖像可以當(dāng)作,由于成像原理不同而造成的二者相同信息的缺失可以看作噪聲干擾?;バ畔⒁话阌糜诿枋鰞蓚€(gè)系統(tǒng)間,即和的統(tǒng)計(jì)相關(guān)性,用熵來(lái)描述為:
式中:()和()分別是系統(tǒng)和的熵;(,)是二者的聯(lián)合熵;(|)為已知系統(tǒng)時(shí)的條件熵;(|)為已知系統(tǒng)時(shí)的條件熵;P()和P()分別是系統(tǒng)和完全獨(dú)立時(shí)的概率分布;P(,)是系統(tǒng)和的聯(lián)合概率分布;P|B(|)是已知系統(tǒng)時(shí)的條件概率分布;P|A(|)是已知系統(tǒng)時(shí)的條件概率分布。
將式(2)、(3)和(4)代入式(1)中,可得兩幅圖像互信息的表達(dá)式:
即在搜索空間中進(jìn)行搜索,當(dāng)互信息達(dá)到最大時(shí),此時(shí)對(duì)應(yīng)的搜索空間的空間位置即為配準(zhǔn)位置。
由于粒子群算法具有實(shí)現(xiàn)容易、精度高、收斂快等優(yōu)點(diǎn),且可以并行計(jì)算,適合在GPU上加速運(yùn)算,因此本文搜索算法采用粒子群算法[9]。
粒子群優(yōu)化算法可表述為:搜索空間即為由方向位移、方向位移、旋轉(zhuǎn)角度、縮放倍數(shù)形成的4維空間,這個(gè)空間中的每個(gè)點(diǎn)為一個(gè)潛在解,被稱(chēng)為“粒子”。每個(gè)粒子有一個(gè)由目標(biāo)函數(shù)確定的適應(yīng)值,適應(yīng)值的大小等于配準(zhǔn)圖像經(jīng)此粒子的變換參數(shù)變換后與待配準(zhǔn)圖像的歸一化互信息的值。每個(gè)粒子的變化量稱(chēng)為速度,一般由如下參數(shù)確定,即此粒子自身的歷史最優(yōu)值、全局變量的歷史最優(yōu)值、慣性權(quán)重、學(xué)習(xí)因子,以及若干隨機(jī)變量,如式(8)所示:
=×+1×rand×(best-pop)+
2×rand (best-pop) (8)
式中:是粒子速度;是慣性權(quán)重;1和2為學(xué)習(xí)因子,一般取1=2=2;best為個(gè)體極值;best為群體極值;rand為0~1之間的一個(gè)隨機(jī)數(shù),pop為某個(gè)粒子所在解空間的位置。慣性權(quán)重較小時(shí),偏向于發(fā)揮粒子的局部搜索能力,反之則偏向于發(fā)揮粒子的全局搜索能力,慣性權(quán)重一般取0.8。通過(guò)多次迭代,粒子在解空間中逐漸收斂到最優(yōu)解,這個(gè)解即為配準(zhǔn)所需的變換參數(shù)。
配準(zhǔn)過(guò)程的流程框圖如圖2所示。
圖2 基于互信息和粒子群算法的流程圖
在這一過(guò)程中,初始化階段包括原始圖像的讀取、灰度化、對(duì)輸入圖像進(jìn)行兩次采樣,以及將Host端的圖像數(shù)據(jù)傳遞到Device端。其中,進(jìn)行的兩次采樣分別生成原圖1/4和1/16的圖像,目的是為了減少運(yùn)算量,即在搜索階段首先在1/16的圖像上進(jìn)行第一階段搜索;然后以輸出參數(shù)作為下一階段搜索的初始狀態(tài)在1/4大小的圖像上進(jìn)行第二階段搜索;最后以第二次搜索結(jié)果為初始狀態(tài),在原圖上進(jìn)行搜索。這一過(guò)程可概括為由粗到精的搜索配準(zhǔn)。圖2所示的過(guò)程為順序執(zhí)行,即邏輯上由Host端控制。
對(duì)于粒子群搜索算法,由于每個(gè)粒子在搜索空間中獨(dú)立計(jì)算,因此并行性主要體現(xiàn)在每一輪迭代時(shí),可以對(duì)個(gè)粒子并行地計(jì)算粒子適應(yīng)度、統(tǒng)計(jì)個(gè)體最優(yōu)適應(yīng)度、粒子速度更新和搜索空間中的位置更新,之后進(jìn)行下一輪迭代。對(duì)于粒子適應(yīng)度的計(jì)算,串行和并行結(jié)構(gòu)的對(duì)比圖如圖3所示。
根據(jù)圖3所示,豎實(shí)線左側(cè)為串行粒子群算法,右側(cè)為相應(yīng)的并行化改寫(xiě),對(duì)應(yīng)部分由橫虛線分割。在并行算法中,每輪迭代個(gè)粒子同時(shí)進(jìn)行運(yùn)算,運(yùn)算完成后,需要柵欄同步函數(shù)進(jìn)行同步,之后才能進(jìn)行下一步計(jì)算。在Device端的顯存中,粒子所對(duì)應(yīng)的圖像數(shù)據(jù)組織成三維結(jié)構(gòu),即每個(gè)粒子對(duì)應(yīng)的圖像保存為二維矩陣,個(gè)粒子對(duì)應(yīng)的圖像數(shù)據(jù)按軸排列為三維數(shù)據(jù)。
對(duì)于每個(gè)粒子,在計(jì)算互信息時(shí)也可進(jìn)行并行優(yōu)化,整個(gè)求取互信息的過(guò)程都充分運(yùn)用GPU的并行執(zhí)行能力,對(duì)于每個(gè)粒子對(duì)應(yīng)圖像進(jìn)行的仿射變換,利用GPU進(jìn)行并行化改進(jìn)。即個(gè)粒子對(duì)應(yīng)的圖像按照各自的變換參數(shù)并行變換,且二次插值時(shí),每個(gè)像素點(diǎn)并行地進(jìn)行插值運(yùn)算。
本實(shí)驗(yàn)的硬件平臺(tái)為:
1)Host端CPU:Intel(R)Core(TM) i3-3220 CPU @3.30GHz;
2)Host端內(nèi)存:4GB;
3)Device端設(shè)備:Nvidia GeForce GTX 650 GPU,時(shí)鐘頻率1072MHz,384個(gè)流處理器,1GB顯存。
軟件平臺(tái)為:
1)編程環(huán)境:VS2010;
2)并行計(jì)算開(kāi)發(fā)工具:CUDA 5.5;
3)圖像處理相關(guān)函數(shù)庫(kù):Opencv 2.4.8;
4)操作系統(tǒng):Win7 64位系統(tǒng)。
對(duì)本文提出的算法進(jìn)行了驗(yàn)證。實(shí)驗(yàn)圖像共3組,每組包含一幅可見(jiàn)光圖像(左)和一幅紅外圖像(右),各組圖像分別如圖4、圖5、圖6所示。
上述3組圖像,尺寸依次為400×400、512×512和639×480。對(duì)于每組圖像,分別將紅外圖像進(jìn)行單獨(dú)水平位移10像素、單獨(dú)垂直位移10像素、單獨(dú)旋轉(zhuǎn)10°、單獨(dú)縮放0.8倍、單獨(dú)縮放1.2倍,以及綜合進(jìn)行位移、旋轉(zhuǎn)和縮放(0.8倍和1.2倍),即每組圖像分別進(jìn)行7組實(shí)驗(yàn),驗(yàn)證算法的有效性。實(shí)驗(yàn)結(jié)果如表1~表3所示。
通過(guò)實(shí)驗(yàn)可知,基于互信息和粒子群的算法可以較為準(zhǔn)確地搜索到可見(jiàn)光和紅外圖像的變換參數(shù)。如實(shí)驗(yàn)數(shù)據(jù)可知,旋轉(zhuǎn)誤差在1°以內(nèi),縮放倍數(shù)誤差在0.05以下,在水平和垂直方向上,誤差為3個(gè)像素以內(nèi),配準(zhǔn)效果較好。以上述試驗(yàn)中的第6組實(shí)驗(yàn)為例,配準(zhǔn)效果如圖7、圖8和圖9所示。
圖3 粒子群串/并行算法對(duì)比圖
Fig. 3 Comparison of particle swarm/parallel algorithm
圖4 機(jī)場(chǎng)A場(chǎng)景可見(jiàn)光(左)與紅外(右)圖像對(duì)比圖
圖5 機(jī)場(chǎng)B場(chǎng)景可見(jiàn)光(左)與紅外(右)圖像對(duì)比圖
圖6 房屋場(chǎng)景可見(jiàn)光(左)與紅外(右)圖像對(duì)比圖
表1 機(jī)場(chǎng)A圖像變換參數(shù)搜索結(jié)果
表2 機(jī)場(chǎng)B圖像變換參數(shù)搜索結(jié)果
表3 房屋圖像變換參數(shù)搜索結(jié)果
如圖7、圖8和圖9所示,3組圖像從左向右依次為參考圖像、待配準(zhǔn)圖像、待配準(zhǔn)圖像變換后的圖像,和將二者進(jìn)行加權(quán)融合的圖像。從最后一幅加權(quán)融合圖像的邊緣部分可以較為清楚地看到配準(zhǔn)結(jié)果較為準(zhǔn)確,結(jié)合表中所示的變換參數(shù),驗(yàn)證算法的正確性。
本組實(shí)驗(yàn)驗(yàn)證了異構(gòu)并行計(jì)算的加速比(加速比=串行運(yùn)算耗時(shí)/異構(gòu)并行運(yùn)算耗時(shí))。對(duì)于本文算法,影響運(yùn)算速度的因素包括圖像大小、粒子規(guī)模,和迭代次數(shù),而與圖像內(nèi)容無(wú)關(guān)。本實(shí)驗(yàn)以機(jī)場(chǎng)A圖像為例,分別單獨(dú)改變其中的某一個(gè)因素,統(tǒng)計(jì)隨這一因素變化情況下加速比的變化情況:
1)圖像尺寸對(duì)加速比的影響根據(jù)算法原理,圖像尺寸的變化會(huì)直接影響仿射變換的執(zhí)行速度,并會(huì)一定程度上影響互信息的統(tǒng)計(jì)速度。由于kernel函數(shù)首次調(diào)用時(shí)可能會(huì)有延時(shí),為了結(jié)果更為準(zhǔn)確,在實(shí)驗(yàn)中對(duì)每幅圖像重復(fù)試驗(yàn)1000次取加速比的平均值。粒子規(guī)模為20,每輪迭代50次,通過(guò)逐漸增大圖像尺寸觀察二者以及配準(zhǔn)算法加速比的變化情況,如圖10、圖11,和圖12所示。
圖10、圖11,和圖12中,橫坐標(biāo)表示圖像邊長(zhǎng),以像素?cái)?shù)為單位,范圍從100×100~10000×10000;縱坐標(biāo)表示加速比。如圖10和圖11所示,仿射變換的加速比可達(dá)10~11倍,互信息的加速比可達(dá)5倍。其原因在于隨著圖像像素?cái)?shù)的增加,仿射變換過(guò)程中的插值次數(shù)也隨之增加,插值算法可采用并行運(yùn)算,即每個(gè)像素的插值過(guò)程相互之間沒(méi)有影響。因此,加速比可以達(dá)到一個(gè)較高的倍數(shù),之后由于硬件資源的限制,加速比趨于穩(wěn)定;而互信息的統(tǒng)計(jì)大體分為兩步,第一步根據(jù)兩幅輸入圖像生成一個(gè)256×256的聯(lián)合灰度直方圖矩陣,第二步由此矩陣統(tǒng)計(jì)得到互信息。由此可知,圖像尺寸的增大主要影響互信息統(tǒng)計(jì)的第一步,對(duì)第二步的影響并不大。此外,第二步通過(guò)歸約求和的算法進(jìn)行并行加速,加速比僅受聯(lián)合灰度直方圖矩陣尺寸的影響,而此矩陣的尺寸為定值。綜上所述,隨圖像尺寸的增大,仿射變換和統(tǒng)計(jì)互信息的加速比都會(huì)增大并趨于平穩(wěn),仿射變換可以達(dá)到更大的加速比。由于仿射變換和統(tǒng)計(jì)互信息過(guò)程中加速比的限制,配準(zhǔn)的整體過(guò)程加速比最大可達(dá)8倍,如圖13所示。
2)粒子規(guī)模對(duì)加速比的影響
考慮到實(shí)驗(yàn)平臺(tái)的內(nèi)存限制,實(shí)驗(yàn)選擇尺寸為400×400的圖像,每輪迭代50次,粒子規(guī)模從10逐漸增大,統(tǒng)計(jì)得到算法加速比如圖14所示。
圖13中,橫坐標(biāo)表示粒子規(guī)模大小??v坐標(biāo)表示對(duì)應(yīng)規(guī)模下配準(zhǔn)算法的加速比。由圖可知,隨著粒子規(guī)模的增大,加速比隨之增大并趨于定值。此例中,加速比最大約為8.5。
3)迭代次數(shù)對(duì)加速比的影響
同樣考慮實(shí)驗(yàn)效率,選擇實(shí)驗(yàn)圖像大小為200×200,粒子規(guī)模為50,逐漸增大迭代次數(shù),統(tǒng)計(jì)配準(zhǔn)算法加速比同迭代次數(shù)的關(guān)系,如圖13所示。圖13中,橫坐標(biāo)為迭代次數(shù),縱坐標(biāo)為相應(yīng)迭代次數(shù)下配準(zhǔn)算法的加速比。由圖可知,在圖像尺寸一定、粒子規(guī)模一定的條件下,配準(zhǔn)算法的加速比隨迭代次數(shù)的增加而增大,最終趨于定值。此例中,加速比最大約為5.5。
圖7 機(jī)場(chǎng)A圖像配準(zhǔn)
Fig.7 Image registration of Airport A
圖8 機(jī)場(chǎng)B圖像配準(zhǔn)
圖9 房屋圖像配準(zhǔn)
圖10 不同尺寸圖像仿射變換的加速比
圖11 不同尺寸統(tǒng)計(jì)兩幅圖像互信息的加速比
Fig.11 The speedup of two images different size statistics mutual information
圖12 不同尺寸配準(zhǔn)算法加速比統(tǒng)計(jì)圖
Fig. 12 The speedup of different size registration algorithms
圖13 不同粒子規(guī)模配準(zhǔn)算法加速比統(tǒng)計(jì)圖
圖14 不同迭代次數(shù)配準(zhǔn)算法加速比統(tǒng)計(jì)圖
本文利用粒子群搜索和互信息相結(jié)合的算法,由粗到精地對(duì)異源圖像進(jìn)行了配準(zhǔn),并利用CUDA C語(yǔ)言對(duì)算法進(jìn)行了異構(gòu)實(shí)現(xiàn),通過(guò)大量實(shí)驗(yàn)驗(yàn)證了算法的魯棒性。對(duì)于異源可見(jiàn)光與紅外圖像配準(zhǔn),以互信息為相似性度量條件,以粒子群算法為搜索算法,在搜索空間中搜索極大值點(diǎn),通過(guò)改變圖像分辨率,由粗到精,逐步實(shí)現(xiàn)可見(jiàn)光圖像與紅外圖像的配準(zhǔn)。實(shí)驗(yàn)證明,在不降低精度的前提下提高了算法效率。針對(duì)粒子群算法、統(tǒng)計(jì)互信息和仿射變換3個(gè)部分的內(nèi)在并行性,利用CUDA C語(yǔ)言實(shí)現(xiàn)GPU-CPU異構(gòu)編程。實(shí)驗(yàn)證明得到了很好的加速比。
[1] 倪國(guó)強(qiáng), 劉瓊. 多源圖像配準(zhǔn)技術(shù)分析與展望[J]. 光電工程, 2004, 31(9): 1-6.
NI Guoqiang, LIU Qiong. Analysis and prospect of multi-source image registration technology[J]., 2004, 31(9): 1-6.
[2] 趙德利, 朱尤攀, 吳誠(chéng), 等. 一種改進(jìn)的聯(lián)合點(diǎn)特征與灰度特征的紅外圖像配準(zhǔn)算法研究[J].紅外技術(shù), 2014, 36(10): 820-826.
ZHAO Deli, ZHU Youpan, WU Chen, et al. An improved algorithm for infrared image registration based on the characteristics of combined points and gray level [J]., 2014, 36 (10): 820-826.
[3] 羅藝, 杜海文, 軒永波, 等. 基于極大似然法的異類(lèi)傳感器配準(zhǔn)方法研究[J]. 電光與控制, 2009, 16(5): 78-80.
LUO Yi, DU Haiwen, XUAN Yongbo, et al. Spatial regeistration of different sensors using maximum likelihood method[J]., 2009, 16 (5): 78-80.
[4] 鄒剛. 基于GPU的圖像配準(zhǔn)及拼接技術(shù)研究[D]. 天津: 天津大學(xué), 2012.
ZOU Gang. Research on image registration and mosaic technology based on GPU [D]. Tianjin: Tianjin University, 2012.
[5] Brown Gottesfeld L. A survey of image registration techniques[J]., 1992, 24(4): 325-376 .
[6] Alahi A, Ortiz R, Vandergheynst P. Freak: fast retina keypoint[C]//2012(CVPR), 2012: 510-517.
[7] 科克, 胡文美, 陳曙暉, 等. 大規(guī)模并行處理器編程實(shí)戰(zhàn)[M]. 北京: 清華大學(xué)出版社, 2010: 1-6.
Kirk, HU Wenmei, CHEN Shuhui, et al.[M]. Beijing: Tsinghua University press, 2010.
[8] Viola P, Wells W M. Alignment by maximization of mutual information[J]., 1997, 24(2): 137-154 .
[9] 陳顯毅. 圖像配準(zhǔn)技術(shù)及其MATLAB編程實(shí)現(xiàn)[M]. 北京: 電子工業(yè)出版社, 2009.
CHEN Xianyi.[M]. Beijing: Publishing House of Electronics Industry, 2009.
Heterogeneous Images Registration Based on Mutual Information and ParticleSwarm Optimization Algorithm Using GPU Parallel Architecture
YU Chunchao1,YANG Zhixiong1,XIA Zongze2,YUAN Xiaochun1,YAN Min1
(1.Kunming Institute of Physics, Kunming 650223, China; 2.State Grid Liaoyang Electric Power Supply Company, Liaoyang 111000, China)
For visible and infrared images registration, in order to get results gradually by changing resolution from low to high, mutual information is used as similarity measure condition, and the particle swarm optimization algorithm is used to search for the maxima in searching space. Experiments show that the algorithm improves the computational efficiency without reducing the accuracy. For the parallelism in particle swarm optimization algorithm, mutual information algorithm, and affine transformation, a heterogeneous GPU-CPU programming is achieved with CUDA C. Experiments show that the algorithm improves the computational efficiency without reducing the accuracy and affects the extreme values of speedup ratio. It proves that the algorithm has a high correct matching rate, high robustness, and high computational efficiency.
Heterogeneous images,mutual information,particle swarm optimization,images registration,GPU
TP391.4
A
1001-8891(2016)06-0938-09
2016-06-21;
2016-08-29.
余春超(1977-),碩士,高級(jí)工程師,主要從事光譜技術(shù)及紅外圖像處理技術(shù),E-mail:beyond_ycc@sina.com。