基于相鄰基準子圖灰度統計相關的快速互信息匹配算法*

符艷軍，孫開鋒

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

0 引言

在基于視覺的飛行器導航與制導、遙感衛星災害監控與環境監測、醫學圖像分析等應用中， 常常需要對不同成像傳感器獲取的異源圖像進行匹配。

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

文中利用匹配過程中相鄰基準子圖間灰度統計特征之間的相關性，通過差量法減少每一個匹配位置互信息的計算量來加快匹配速度。與已有的以灰度壓縮、特征提取或采取優化搜索策略等快速互信息匹配方法相比，該方法既沒有減少參與匹配的像素數，也無需對圖像灰度等級進行處理，不會對匹配精度造成影響。

1 圖像熵及互信息匹配

圖像熵描述了圖像信源的平均信息量，反映了圖像灰度的統計信息，相似的兩幅圖像其圖像熵也相近。利用互信息進行圖像匹配的實質是：當兩幅圖像在空間位置配準時，其重疊部分所對應像素對的灰度互信息達到最大值。 兩幅圖像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)表示兩幅圖像的聯合熵。熵通常由變量的概率密度來表示。針對灰度圖像，可以分別用兩幅圖像的灰度直方圖和聯合直方圖進行估計。 計算公式分別為：

(2)

(3)

(4)

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

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

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

(5)

在飛行器導航中，匹配的目的是確定飛行器的當前位置，通過計算實測圖與基準圖中每一個匹配位置對應基準子圖的互信息并找出最大互信息值對應位置，即可確定基準圖與實測圖的相對位置, 實現對飛行器的定位。

2 快速互信息匹配算法

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

圖1 相鄰基準子圖及其灰度直方圖

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

1)獲取實測圖像A和基準圖像R，將兩幅圖像的灰度值調整到同一灰度區間。

圖2 文中匹配算法流程圖

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

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

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

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

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

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

(6)

(7)

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

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

(8)

④計算實測圖與當前基準子圖的聯合直方圖矩陣hAS，找出hAS和前一基準子圖聯合直方圖矩陣hAS0的不同元素，分別組成矩陣hΔAS和hΔAS0，按照式(8)、式(9)分別計算hΔAS和hΔAS0對應的聯合熵ΔAS及ΔAS0；其中，normAS0為與hΔAS0中元素對應的前一基準子圖聯合熵計算過程中的求和項，因為這些求和項在前一匹配位置已經計算過，此次無需再重復計算，而只需計算這些項的和即可。

(9)

(10)

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

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

(11)

⑥按照式(12)計算實測圖與當前基準子圖的歸一化互信息。

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

(12)

6) 遍歷整個基準圖，選取歸一化互信息最大值對應的位置作為最終匹配點。

3 算法時間復雜性分析

文中算法與傳統方法相比，計算量的差別主要體現在每一個匹配位置處基準子圖各統計量的計算上，為此只分析基準子圖各統計量計算的時間復雜度。

設基準圖尺寸為N×N，實測圖尺寸為n×n，則基準子圖大小與實測圖大小一致均為n×n。

1)基準子圖直方圖計算

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

2)基準子圖熵及聯合熵計算

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

為此，文中在標準圖像集上進行了大量的仿真實驗，通過對仿真結果進行統計得出：相鄰子圖間灰度直方圖中發生變化的元素平均為50%左右，聯合灰度直方圖中發生變化的元素平均為75%左右，這就意味著文中方法相比傳統方法，其基準子圖熵值的計算量減少了近50%，聯合熵計算量減少了約25%，時間復雜度或時間頻度有所降低。部分實驗結果列于表1～表4。

另一方面，因文中方法需要存儲前一匹配位置基準子圖的多個統計量，算法的空間復雜度有所增加，但由于硬件技術的發展大大提高了計算機的存儲容量，使得存儲容量的局限性對于算法的影響大大降低。

4 仿真實驗

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

匹配過程中設: 基準子圖直方圖計算中需要統計的像素灰度值的次數為M，基準子圖熵值計算中需要計算的求和項的項數為M1，基準子圖與實測圖聯合熵計算中需要計算的求和項項數為M2。

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

圖3 SAR與可見光圖像匹配

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

5 結語

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

表1 SAR與可見光圖像在列方向匹配時部分位置計算量統計結果

表2 SAR與可見光圖像在行方向匹配時部分位置計算量統計結果

表3 差量法與傳統標準法在SAR與可見光圖像上的匹配結果

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