基于相鄰基準(zhǔn)子圖灰度統(tǒng)計(jì)相關(guān)的快速互信息匹配算法*

符艷軍，孫開鋒

(1 西安外事學(xué)院，西安 710077； 2 西安精密機(jī)械研究所，西安 710075)

0 引言

在基于視覺的飛行器導(dǎo)航與制導(dǎo)、遙感衛(wèi)星災(zāi)害監(jiān)控與環(huán)境監(jiān)測、醫(yī)學(xué)圖像分析等應(yīng)用中， 常常需要對不同成像傳感器獲取的異源圖像進(jìn)行匹配。

基于互信息測度的匹配方法是一種完全基于圖像灰度統(tǒng)計(jì)概率的匹配方法，不需要對原圖像間的灰度關(guān)系作任何假設(shè)，非常適合于異源圖像的匹配，但互信息測度最大的缺陷是計(jì)算量大， 匹配耗時(shí)較長[1]。為此，研究者們就圍繞減少互信息測度計(jì)算量以及改進(jìn)搜索策略等方面提出了相應(yīng)的加速方法。文獻(xiàn)[2-5]分別采用模擬退火、粒子群等優(yōu)化搜索算法提高匹配速度；文獻(xiàn)[6]采用灰度壓縮和粒子群優(yōu)化算法相結(jié)合加快匹配速度；文獻(xiàn)[7]通過提取多層次特征點(diǎn)減少了互信息測度的計(jì)算量；文獻(xiàn)[8-10]采用由粗到精的多層匹配方法對算法進(jìn)行加速；文獻(xiàn)[11-15]將分層匹配與各種優(yōu)化算法相結(jié)合提高匹配速度。

文中利用匹配過程中相鄰基準(zhǔn)子圖間灰度統(tǒng)計(jì)特征之間的相關(guān)性，通過差量法減少每一個(gè)匹配位置互信息的計(jì)算量來加快匹配速度。與已有的以灰度壓縮、特征提取或采取優(yōu)化搜索策略等快速互信息匹配方法相比，該方法既沒有減少參與匹配的像素?cái)?shù)，也無需對圖像灰度等級進(jìn)行處理，不會對匹配精度造成影響。

1 圖像熵及互信息匹配

圖像熵描述了圖像信源的平均信息量，反映了圖像灰度的統(tǒng)計(jì)信息，相似的兩幅圖像其圖像熵也相近。利用互信息進(jìn)行圖像匹配的實(shí)質(zhì)是：當(dāng)兩幅圖像在空間位置配準(zhǔn)時(shí)，其重疊部分所對應(yīng)像素對的灰度互信息達(dá)到最大值。 兩幅圖像U和V的互信息I(U,V)可以用下式表示：

I(U,V)=H(U)+H(V)-H(U,V)

(1)

其中:H(U)、H(V) 分別表示圖像U及V的熵；H(U,V)表示兩幅圖像的聯(lián)合熵。熵通常由變量的概率密度來表示。針對灰度圖像，可以分別用兩幅圖像的灰度直方圖和聯(lián)合直方圖進(jìn)行估計(jì)。 計(jì)算公式分別為：

(2)

(3)

(4)

式(2)～式(4)中，hU(i)表示圖像U中灰度值為i的像素出現(xiàn)的次數(shù)；hV(i)表示圖像V中灰度值為i的像素出現(xiàn)的次數(shù)；hUV(i,j)表示同一像素位置上圖像U中灰度值為i、圖像V中灰度值為j的像素對出現(xiàn)的次數(shù)；T為圖像的總像素?cái)?shù)。

研究發(fā)現(xiàn)[16]，互信息本身的大小與待配準(zhǔn)兩圖像間的重疊度有一定的關(guān)聯(lián)性，為了消除這種關(guān)聯(lián)關(guān)系，文獻(xiàn)[17]提出了利用歸一化互信息(normalized mutual information， NMI)作為相似性測度。NMI能夠減少對圖像重疊部分的敏感性，實(shí)驗(yàn)表明它比標(biāo)準(zhǔn)的互信息方法更具魯棒性。NMI表達(dá)式為：

NMI(U,V)=(H(U)+H(V))/H(U,V)

(5)

在飛行器導(dǎo)航中，匹配的目的是確定飛行器的當(dāng)前位置，通過計(jì)算實(shí)測圖與基準(zhǔn)圖中每一個(gè)匹配位置對應(yīng)基準(zhǔn)子圖的互信息并找出最大互信息值對應(yīng)位置，即可確定基準(zhǔn)圖與實(shí)測圖的相對位置, 實(shí)現(xiàn)對飛行器的定位。

2 快速互信息匹配算法

與基于特征的圖像匹配算法相比，基于最大歸一化互信息的圖像匹配最大缺點(diǎn)是計(jì)算速度慢。 每計(jì)算一次互信息測度，需要遍歷圖像中的所有像素，這使得計(jì)算量大大增加。由式(5)可以看出，歸一化互信息的計(jì)算量主要集中在兩幅圖像各自的熵及其聯(lián)合熵的計(jì)算。從式(2)～式(4)可以看出，圖像熵及聯(lián)合熵均為求和運(yùn)算，其計(jì)算量主要集中于每一個(gè)求和項(xiàng)的計(jì)算，因此，如果能減少每次匹配過程中需要計(jì)算的求和項(xiàng)的項(xiàng)數(shù)， 則可以加快匹配速度。因?yàn)閳D像分布具有塊狀結(jié)構(gòu)，在遍歷式搜索匹配過程中，各相鄰基準(zhǔn)子圖的灰度統(tǒng)計(jì)量具有很強(qiáng)的相關(guān)性，基準(zhǔn)圖像(i,j)位置處對應(yīng)子圖的灰度直方圖與其四鄰域位置處對應(yīng)子圖的灰度直方圖非常相近。

圖1 相鄰基準(zhǔn)子圖及其灰度直方圖

如圖1所示，圖1(a)為原始基準(zhǔn)圖像，圖1(b)為圖1(a)的一個(gè)子圖，圖1(c)為圖1(b)在原基準(zhǔn)圖上右移一個(gè)像素對應(yīng)的子圖，圖1(d)和圖1(e)分別為圖1(b)和圖1(c)的灰度直方圖。從圖中可以看出，相鄰兩個(gè)子圖的直方圖非常相似，也即兩相鄰子圖存在大量的灰度值其出現(xiàn)的頻率沒有發(fā)生變化，相應(yīng)的，在式(2)～式(4)的求和項(xiàng)中分別有很多項(xiàng)數(shù)沒有發(fā)生變化，因此， 當(dāng)前匹配位置處基準(zhǔn)子圖的熵以及與實(shí)測圖的聯(lián)合熵的計(jì)算可以前一匹配位置處基準(zhǔn)子圖的熵及與實(shí)測圖的聯(lián)合熵的計(jì)算為基準(zhǔn)，通過比較相鄰兩個(gè)匹配位置處基準(zhǔn)子圖對應(yīng)的直方圖及聯(lián)合直方圖的變化情況，用差量法進(jìn)行計(jì)算。即在整個(gè)匹配過程中，除了要完整計(jì)算基準(zhǔn)圖上第一匹配位置(1,1)處對應(yīng)基準(zhǔn)子圖的熵以及與實(shí)測圖的聯(lián)合熵外，其它所有匹配位置的基準(zhǔn)子圖的熵及與實(shí)測圖的聯(lián)合熵的計(jì)算主要集中在與前一匹配位置相比灰度直方圖矩陣和聯(lián)合直方圖矩陣中發(fā)生變化的那些元素對應(yīng)的求和項(xiàng)上。據(jù)此， 文中互信息匹配算法的流程圖如圖2所示。具體步驟如下:

1)獲取實(shí)測圖像A和基準(zhǔn)圖像R，將兩幅圖像的灰度值調(diào)整到同一灰度區(qū)間。

圖2 文中匹配算法流程圖

2)按照式(2)或式(3)計(jì)算A的信息熵H(A)并存儲。

3)從基準(zhǔn)圖像R的左上角(1,1)點(diǎn)截取與實(shí)測圖像A大小相等的第1幅基準(zhǔn)子圖S1,統(tǒng)計(jì)S1各灰度值出現(xiàn)次數(shù)并存入初始灰度直方圖矩陣hS1，按照式(2)或式(3)計(jì)算S1的信息熵H(S1)并存儲，最后將hS1和H(S1)分別賦值給hS0和H(S0)。

4)統(tǒng)計(jì)實(shí)測圖A與第1幅基準(zhǔn)子圖S1對應(yīng)像素位置的灰度對出現(xiàn)次數(shù)并存入初始聯(lián)合直方圖矩陣hAS1，按照式(4)計(jì)算初始聯(lián)合熵H(AS1)并保存，將hAS1和H(AS1)分別賦值給hAS0和H(AS0)；再按照式(5)計(jì)算實(shí)測圖與基準(zhǔn)子圖的歸一化互信息值并存儲。

5)從基準(zhǔn)圖上(1,2)點(diǎn)位置開始，逐行逐像素按照“Z”字形進(jìn)行迭代匹配。匹配過程中， 當(dāng)前匹配位置記為(i,j)，當(dāng)j>1時(shí)，當(dāng)前匹配位置的前一匹配位置應(yīng)為(i,j-1)；當(dāng)j=1時(shí)，則當(dāng)前匹配位置的前一匹配位置應(yīng)為(i-1,j)，則當(dāng)前位置(i,j)處的歸一化互信息可按如下步驟計(jì)算:

①通過比較當(dāng)前匹配位置(i,j)處基準(zhǔn)子圖和前一匹配位置(i,j-1)或(i-1,j)處基準(zhǔn)子圖對應(yīng)列或行像素變化情況，按照差量法將前一基準(zhǔn)子圖的直方圖矩陣hS0更新為當(dāng)前子圖的直方圖矩陣hS。具體做法為: 在行方向匹配時(shí)，當(dāng)前匹配位置的前一匹配位置為(i,j-1)，則以hS0為基準(zhǔn)，減去前一基準(zhǔn)子圖第一列對應(yīng)像素灰度值出現(xiàn)次數(shù)，加上當(dāng)前基準(zhǔn)子圖最后一列對應(yīng)像素灰度值出現(xiàn)次數(shù)，即得到hS；在列方向匹配時(shí)，當(dāng)前匹配位置的前一匹配位置為(i-1,j)，則以hS0為基準(zhǔn)，減去前一基準(zhǔn)子圖第一行對應(yīng)像素灰度值出現(xiàn)次數(shù)，加上當(dāng)前基準(zhǔn)子圖最后一行對應(yīng)像素灰度值出現(xiàn)次數(shù)，即得到hS。

②找出hS與hS0的不同元素，分別組成矩陣hΔS和hΔS0，按照式(6)、式(7)分別計(jì)算hΔS和hΔS0對應(yīng)的信息熵ΔS及ΔS0；其中，normS0(S)為與hΔS0中元素對應(yīng)的前一基準(zhǔn)子圖熵值計(jì)算過程中的求和項(xiàng)，因?yàn)檫@些求和項(xiàng)在前一匹配位置已經(jīng)計(jì)算過，此次無需再重復(fù)計(jì)算，而只需計(jì)算這些項(xiàng)的和即可。

(6)

(7)

③以前一基準(zhǔn)子圖的信息熵H(S0)為基準(zhǔn)，按照式(8)通過差量法將H(S0)更新為當(dāng)前基準(zhǔn)子圖的信息熵H(S)。

H(S)=H(S0)-ΔS0+ΔS

(8)

④計(jì)算實(shí)測圖與當(dāng)前基準(zhǔn)子圖的聯(lián)合直方圖矩陣hAS，找出hAS和前一基準(zhǔn)子圖聯(lián)合直方圖矩陣hAS0的不同元素，分別組成矩陣hΔAS和hΔAS0，按照式(8)、式(9)分別計(jì)算hΔAS和hΔAS0對應(yīng)的聯(lián)合熵ΔAS及ΔAS0；其中，normAS0為與hΔAS0中元素對應(yīng)的前一基準(zhǔn)子圖聯(lián)合熵計(jì)算過程中的求和項(xiàng)，因?yàn)檫@些求和項(xiàng)在前一匹配位置已經(jīng)計(jì)算過，此次無需再重復(fù)計(jì)算，而只需計(jì)算這些項(xiàng)的和即可。

(9)

(10)

⑤以前一基準(zhǔn)子圖的聯(lián)合熵H(AS0)為基準(zhǔn)，按照式(11)通過差量法將H(AS0)更新為當(dāng)前基準(zhǔn)子圖的聯(lián)合熵H(AS)。

H(AS)=H(AS0)-ΔAS0+ΔAS

(11)

⑥按照式(12)計(jì)算實(shí)測圖與當(dāng)前基準(zhǔn)子圖的歸一化互信息。

NMI(A,S)=(H(A)+H(S))/H(A,S)

(12)

6) 遍歷整個(gè)基準(zhǔn)圖，選取歸一化互信息最大值對應(yīng)的位置作為最終匹配點(diǎn)。

3 算法時(shí)間復(fù)雜性分析

文中算法與傳統(tǒng)方法相比，計(jì)算量的差別主要體現(xiàn)在每一個(gè)匹配位置處基準(zhǔn)子圖各統(tǒng)計(jì)量的計(jì)算上，為此只分析基準(zhǔn)子圖各統(tǒng)計(jì)量計(jì)算的時(shí)間復(fù)雜度。

設(shè)基準(zhǔn)圖尺寸為N×N，實(shí)測圖尺寸為n×n，則基準(zhǔn)子圖大小與實(shí)測圖大小一致均為n×n。

1)基準(zhǔn)子圖直方圖計(jì)算

每一個(gè)匹配位置處，傳統(tǒng)算法需要遍歷整個(gè)基準(zhǔn)子圖來計(jì)算該子圖的直方圖，而文中方法只須遍歷基準(zhǔn)子圖的第一行(或列)和最后一行(或列)即可獲得該子圖直方圖。對于整個(gè)匹配過程來說，傳統(tǒng)方法其時(shí)間復(fù)雜度為O((N-n)2n2)，文中方法時(shí)間復(fù)雜度為O((N-n)2n)。

2)基準(zhǔn)子圖熵及聯(lián)合熵計(jì)算

由于相鄰基準(zhǔn)子圖其直方圖發(fā)生變化情況跟圖像內(nèi)容有關(guān)，很難從理論上推導(dǎo)兩子圖間有多少個(gè)灰度值其出現(xiàn)頻率會發(fā)生變化，也很難推導(dǎo)出其聯(lián)合直方圖中有多少個(gè)灰度值對的頻率會發(fā)生變化。以直方圖為例，最壞情況時(shí)有2n(n為基準(zhǔn)子圖各邊長所含像素個(gè)數(shù))個(gè)灰度值頻率發(fā)生變化，且2n不超過基準(zhǔn)圖灰度等級數(shù)。

為此，文中在標(biāo)準(zhǔn)圖像集上進(jìn)行了大量的仿真實(shí)驗(yàn)，通過對仿真結(jié)果進(jìn)行統(tǒng)計(jì)得出：相鄰子圖間灰度直方圖中發(fā)生變化的元素平均為50%左右，聯(lián)合灰度直方圖中發(fā)生變化的元素平均為75%左右，這就意味著文中方法相比傳統(tǒng)方法，其基準(zhǔn)子圖熵值的計(jì)算量減少了近50%，聯(lián)合熵計(jì)算量減少了約25%，時(shí)間復(fù)雜度或時(shí)間頻度有所降低。部分實(shí)驗(yàn)結(jié)果列于表1～表4。

另一方面，因文中方法需要存儲前一匹配位置基準(zhǔn)子圖的多個(gè)統(tǒng)計(jì)量，算法的空間復(fù)雜度有所增加，但由于硬件技術(shù)的發(fā)展大大提高了計(jì)算機(jī)的存儲容量，使得存儲容量的局限性對于算法的影響大大降低。

4 仿真實(shí)驗(yàn)

為了驗(yàn)證算法的有效性，將文中方法應(yīng)用于SAR與可見光圖像匹配，以Matlab R2015為仿真平臺，進(jìn)行兩類實(shí)驗(yàn)：一類是通過比較每個(gè)匹配位置處文中方法與傳統(tǒng)遍歷式互信息匹配方法在計(jì)算量方面的差別以及總耗時(shí)來說明文中方法的有效性；另一類是將文中方法與多分辨率分層匹配相結(jié)合，與單純的多分辨率方法進(jìn)行匹配耗時(shí)及匹配精度的比較，其中的多分辨率分層匹配采用d3小波基，小波分解級數(shù)為1。限于篇幅，文中只列出部分結(jié)果。

匹配過程中設(shè): 基準(zhǔn)子圖直方圖計(jì)算中需要統(tǒng)計(jì)的像素灰度值的次數(shù)為M，基準(zhǔn)子圖熵值計(jì)算中需要計(jì)算的求和項(xiàng)的項(xiàng)數(shù)為M1，基準(zhǔn)子圖與實(shí)測圖聯(lián)合熵計(jì)算中需要計(jì)算的求和項(xiàng)項(xiàng)數(shù)為M2。

選取可見光與SAR圖像進(jìn)行匹配，如圖3所示，其中圖3(a)為可見光基準(zhǔn)圖，大小為422×358；圖3(b)為SAR實(shí)測圖像1，大小為150×150，圖3(c)為SAR實(shí)測圖像2，大小為100×90，圖3(d)為基于文中差量法匹配時(shí)實(shí)測圖在基準(zhǔn)圖上的匹配定位結(jié)果；表1為差量法與傳統(tǒng)標(biāo)準(zhǔn)算法列方向匹配時(shí)部分匹配位置計(jì)算量統(tǒng)計(jì)結(jié)果的比較，表2為差量法與傳統(tǒng)標(biāo)準(zhǔn)算法行方向匹配時(shí)部分匹配位置計(jì)算量統(tǒng)計(jì)結(jié)果的比較，表3為文中差量法與傳統(tǒng)標(biāo)準(zhǔn)方法匹配結(jié)果和匹配耗時(shí)的比較，表4為差量法結(jié)合多分辨率分層匹配法與單純多分辨率分層匹配法的匹配結(jié)果和匹配耗時(shí)的比較。

圖3 SAR與可見光圖像匹配

由圖3可以看出，采用歸一化互信息作為匹配測度，可以實(shí)現(xiàn)異源圖像之間的匹配；由表1～表3可以看出，相比傳統(tǒng)方法，文中差量法能夠有效減少除第一個(gè)匹配位置外其它所有匹配位置互信息測度的計(jì)算量，總匹配耗時(shí)平均減少20%以上；由表4可以看出，文中差量法可以與多分辨率分層匹配技術(shù)相結(jié)合，其匹配精度與單純的多分辨率分層匹配算法精度相同(在有些情況下匹配精度有所下降)，但其匹配耗時(shí)平均下降20%左右。

5 結(jié)語

與已有的以灰度壓縮、特征提取或采取優(yōu)化搜索策略等快速互信息匹配方法不同,文中基于圖像的塊狀結(jié)構(gòu)特性以及匹配過程中相鄰子圖統(tǒng)計(jì)特性的相關(guān)性, 通過差量法減少互信息測度的計(jì)算量來加快匹配速度, 其互信息測度計(jì)算精度及匹配精度與原始的基于標(biāo)準(zhǔn)互信息計(jì)算方法的計(jì)算精度及匹配精度相同, 但計(jì)算量大為減少。 另外,如果參與匹配的兩幅圖像為高分辨率清晰圖像，則可以把文中方法與灰度壓縮、多分辨率分層匹配等方法相結(jié)合，即先對待匹配的圖像進(jìn)行灰度壓縮或多分辨率分解，然后采用遍歷式方法進(jìn)行搜索，而在每一個(gè)搜索點(diǎn)處采用文中的差量法計(jì)算互信息測度，這樣可以進(jìn)一步減少匹配耗時(shí)。

表1 SAR與可見光圖像在列方向匹配時(shí)部分位置計(jì)算量統(tǒng)計(jì)結(jié)果

表2 SAR與可見光圖像在行方向匹配時(shí)部分位置計(jì)算量統(tǒng)計(jì)結(jié)果

表3 差量法與傳統(tǒng)標(biāo)準(zhǔn)法在SAR與可見光圖像上的匹配結(jié)果

表4 差量法結(jié)合多分辨率分層算法與單純多分辨率分層算法在SAR與可見光圖像上的匹配結(jié)果