相艷，桂鵬，王碩，許春榮，邵黨國，劉利軍，湯守國

(昆明理工大學(xué) 信息工程與自動化學(xué)院，昆明 650051)

1 引 言

圖像配準(zhǔn)是指通過某種空間變換，能使兩幅圖像的特征點或者對應(yīng)點在空間上達到一致。其中保持不動的圖像稱為參考圖像，進行空間變換的稱為浮動圖像。圖像配準(zhǔn)一般由4個基本模塊組成，即幾何變換、圖像插值、相似性測度和函數(shù)優(yōu)化[1]。相似性測度衡量了參考圖像和浮動圖像匹配的程度。用于圖像配準(zhǔn)的相似性測度有相關(guān)法、最小化灰度差法、灰度比的方差最小化法、劃分一致法、最小化聯(lián)合熵方法和互信息法等，其中歸一化互信息(normalized mutual information，NMI)被認(rèn)為是配準(zhǔn)精度和魯棒性較好的回溯式配準(zhǔn)方法之一[2]。但多模圖像配準(zhǔn)中NMI函數(shù)曲線容易出現(xiàn)局部極值，從而影響配準(zhǔn)的精度和魯棒性。交叉累計剩余熵(cross cumulative residual entropy，CCRE)使用累計剩余分布函數(shù)代替概率密度函數(shù)，能得到更光滑的配準(zhǔn)目標(biāo)函數(shù)[3]，近年來在多模圖像配準(zhǔn)中得到了很好的應(yīng)用。但CCRE 需要計算分布函數(shù)，每次迭代時間較長，影響了配準(zhǔn)效率[4]。Pickering 等人首次提出使用條件方差和(the sum of conditional variance，SCV)進行多模圖像配準(zhǔn),使用高斯牛頓法求解量化后圖像的SCV極小值，從而得到空間變換參數(shù)[5]。文獻[6]的實驗結(jié)果表明這種新的相似性測度與互信息相比，準(zhǔn)確度和魯棒性都有改善。此后，SCV更多被用于視覺追蹤、視覺伺服[7-8]。

SCV的主要缺點是它僅使用參考和浮動圖像的量化信息計算聯(lián)合直方圖[7]。本研究認(rèn)為，在計算SCV的聯(lián)合直方圖時，就應(yīng)該把浮動圖像非網(wǎng)格點的插值方法結(jié)合進去，從而克服SCV的這個限制。浮動圖像非網(wǎng)格點的插值對NMI、CCRE、SCV這些相似性測度曲線的光滑程度都有不同程度的影響。傳統(tǒng)的插值方法有最近鄰插值、線性插值、PV插值，這些插值方法明顯缺陷是會產(chǎn)生插值偽影。插值偽影會產(chǎn)生局部極值，引起全局最優(yōu)值的移動，從而給函數(shù)尋優(yōu)帶來困難，影響配準(zhǔn)結(jié)果[9]。為了減少插值偽影，提高配準(zhǔn)效率，新的插值方法不斷涌現(xiàn)。其中一類方法是改變參與計算聯(lián)合概率密度的鄰域點范圍，并使用改進的核函數(shù)來分配鄰域點的權(quán)重，針對該類方法，Rohde等人通過頻域分析指出，擴大鄰域點范圍、采用高階的插值函數(shù)能從整體上降低局部極值的出現(xiàn)[10]。然而，隨著插值函數(shù)階次的增加，算法的復(fù)雜度也相應(yīng)增加。為此，本研究設(shè)計了一種新的插值函數(shù)來替代高階插值函數(shù)，從而提高現(xiàn)有SCV算法的性能。

實驗結(jié)果表明，本研究方法在進行多模圖像配準(zhǔn)時，在保證準(zhǔn)確度的同時，對空間變換較大和噪聲較強的情況具有較好的魯棒性，比NMI、CCRE和原SCV方法具有更高的配準(zhǔn)成功率。

2 理論與方法2.1 SCV原理

參考圖像和浮動圖像的灰度值可以看作兩個隨機矩陣R、T(n×m)。令x=(x，y)為包含像素坐標(biāo)的向量，(x，y)∈{1，2，…，m}×{1，2，…，n}，w(x，pa)是一個轉(zhuǎn)換函數(shù)，其中空間變換參數(shù)pa將像素x從T(x)映射到R，x→w(x，pa)。參考圖像和浮動圖像之間的條件方差和計算公式為[11]：

(1)

(2)

ε(·)表示數(shù)學(xué)期望運算。

2.2 改進的插值方法

當(dāng)參數(shù)pa將像素x從T(x)映射到R時，該像素經(jīng)過空間變換映射到參考圖像中時通常為非網(wǎng)格點，而非網(wǎng)格點的空間插值方法會影響聯(lián)合直方圖的計算，從而影響SCV曲線的平滑程度，最終影響配準(zhǔn)結(jié)果。常用的PV插值法是根據(jù)權(quán)重分配原則，將每對像素對聯(lián)合直方圖的貢獻按權(quán)重分配給相鄰的四個像素對。PV插值的一個缺陷是當(dāng)插值點正好落在網(wǎng)格點處時，則該網(wǎng)格點處權(quán)重為1，而其他相鄰的三個網(wǎng)格點權(quán)重為0，這會導(dǎo)致相似性測度容易在整數(shù)倍位移處出現(xiàn)局部極值。而改進的PV插值方法通過擴大鄰域點范圍，使用改進的插值函數(shù)來分配鄰域點的權(quán)重。

以2D圖像為例，設(shè)浮動圖像上某點經(jīng)空間變換后與參考圖像的點q(非網(wǎng)格點)對應(yīng)，點q的鄰域點為qi(i=1，2，…，16)，見圖1。p、q點對聯(lián)合直方圖貢獻的概率計算如下：

h(T(p)，R(q))+=ωi，i=1，2，…，16

(3)

其中ωi為每個領(lǐng)域點qi的權(quán)重。

wi=f(dix)·f(diy)

(4)

上式中dix和diy為點qi(i=1，2，…，16)與q在x和y方向上的距離。

圖1 改進的插值方法2D示意圖

文獻[12]指出，插值函數(shù)應(yīng)該滿足以下三個條件：

(1)f(x)≥0；

(3)f(x) 為單調(diào)遞減函數(shù)。

本研究提出了一種滿足上述三個條件的插值函數(shù)，具有如下形式：

(5)

如果令鄰域點與插值點達到最大距離時不分配權(quán)重，即x=2時f(x)=0，則c=2。此時式(5)可簡化為：

(6)

當(dāng)點q落在非網(wǎng)格點上時，參與聯(lián)合直方圖統(tǒng)計的鄰域點數(shù)目有16個；而當(dāng)點q落在網(wǎng)格點q6上時，參與聯(lián)合直方圖統(tǒng)計的鄰域點數(shù)目是最少的，有9 個。能分配到權(quán)重的最多鄰域點數(shù)與最少鄰域點數(shù)的比例為16:9，比傳統(tǒng)PV插值法的比例4:1要小。該比例越小，說明聯(lián)合直方圖的概率分布越分散[12]，從而SCV出現(xiàn)局部極值的幾率降低。

在計算鄰域點權(quán)重時，上式中的公倍數(shù)1/4對鄰域點的影響是相同的，因此可將其提出，公式(6)可以進一步簡化為:

(7)

利用公式(7)計算插值函數(shù)的值，相當(dāng)于在利用公式(6)計算出鄰域點單個方向上的插值權(quán)重之后，再給每一個插值點權(quán)重都擴大4倍。因為這種變換是統(tǒng)一的線性比例變換，因此，不會改變公式(6)的權(quán)重分配，還能提高算法的效率。

2.3 本研究的相似性測度

本研究設(shè)計的相似性測度計算步驟如下：

(1) 對于尋優(yōu)過程中給出的空間變換參數(shù)pa，將像素x從T(x)映射到R。根據(jù)映射至R中的插值點位置計算該點與其16個鄰域點的距離，從而利用公式(7)計算出插值函數(shù)值，進一步利用公式(4)計算出每個鄰域點的權(quán)重。

(2) 根據(jù)公式(3)計算出每一個插值點貢獻給聯(lián)合直方圖的概率，從而得到參考圖像和浮動圖像的聯(lián)合直方圖，用H(u，v)表示，其中u表示參考圖像的灰度，v表示浮動圖像的灰度。u=1，2，…，dR；v=1，2，…dT，dR和dT是R和T的最大像素灰度級數(shù)目。

(3) 將公式(2)中的數(shù)學(xué)期望進一步展開為以下公式進行計算：

(8)

(9)

3 四種相似性測度曲線的比較

以BrainWeb數(shù)據(jù)集[13]中正常腦部磁共振T1加權(quán)第50斷層圖像為參考圖像，見圖2(a)，浮動圖像為相應(yīng)斷層的PD加權(quán)像，見圖2(b)。此外，圖2(c)為根據(jù)公式(8)計算得到的浮動圖像的期望值，圖2(d)為圖1(a)與(c)相減得到的圖像,計算圖2(d)中像素灰度的平方和，即得到(a)與(b)的SCV。

將圖2(b)沿x軸及y軸平移-5～5像素，平移步長為0.1像素，計算其與圖2(a)的NMI、CCRE、SCV以及本研究的相似性測度，得到見圖3結(jié)果(SCV和本研究方法都取負(fù)值)。從圖3中可以看出，NMI、CCRE、SCV在網(wǎng)格點位置都有不同程度的波動，其中NMI的局部極值最為明顯，CCRE次之，SCV也有一定程度的波動，而本研究提出的方法最為平滑。NMI的抖動最為嚴(yán)重，極易使尋優(yōu)陷入局部極值中，導(dǎo)致配準(zhǔn)失敗。本研究方法相比其他三種相似性測度，能保證優(yōu)化算法順利找到全局最值(最大或最小)，從而成功配準(zhǔn)。

圖2 實驗所用的醫(yī)學(xué)圖像

4 多模圖像配準(zhǔn)實驗

本研究分別針對多模圖像配準(zhǔn)中，空間變換范圍不同或噪聲強度不同時，四種相似性測度的表現(xiàn)進行測試，實驗一針對的是空間變換范圍不同的情況，參考圖像為BrainWeb數(shù)據(jù)集中正常腦部磁共振T1加權(quán)第90斷層圖像，浮動圖像為相應(yīng)斷層的T2加權(quán)像并做一定的剛體變換，設(shè)tx、ty分別為浮動圖像沿x和y軸的平移，θ為繞圖像中心點的順時針旋轉(zhuǎn)角度。在一定取值范圍內(nèi)隨機選取三個參數(shù)的值，每一組測試中隨機選取20個不同的剛體變換T=[tx，ty，θ]進行測試，其中第一組測試中tx、ty的絕對值范圍為1～2，θ的絕對值范圍為0～1；第二組中tx、ty的絕對值范圍為14～15，θ的絕對值范圍為2～3；第三組中tx、ty的絕對值范圍加大至19～20，θ的絕對值范圍為3～4,使用步長為5的Powell算法對NMI、CCRE、SCV和本研究方法相似性測度進行尋優(yōu)，得到不同的配準(zhǔn)參數(shù)T′=[tx′，ty′，θ′]，并與真實的參數(shù)T比較，求取與T的誤差△T=[△tx，△ty，△θ]。如果尋優(yōu)陷入局部極值，則配準(zhǔn)失敗。本研究中配準(zhǔn)的性能評估采用以下公式計算：

(10)

如果兩幅圖像完全配準(zhǔn)，那么Err=0，如果Powell尋優(yōu)陷入局部極值，則配準(zhǔn)失敗。圖4為每組測試中20次配準(zhǔn)的成功率。

圖3 NMI、CCRE、SCV和本研究的相似性測度

由上圖可知，第一組測試中，當(dāng)空間變換不大時， NMI、CCRE、SCV和本研究的方法成功率均達到了100%。第二組測試中，當(dāng)空間變換加大時，本研究方法配準(zhǔn)成功率為100%，SCV為89.47%，CCRE為89.47%，NMI為0。第三組測試中，本研究方法配準(zhǔn)成功率為73.68%，SCV為57.89%，CCRE為68.42%，NMI為0。可以看出，當(dāng)浮動圖像相比參考圖像的空間變化增加時，NMI、CCRE、SCV方法配準(zhǔn)成功率都會降低，本研究提出的方法仍然可以達到較高的配準(zhǔn)成功率。

圖4 實驗一的配準(zhǔn)結(jié)果

實驗二針對的是噪聲強度不同的情況。參考和浮動圖像與實驗一相同，但浮動圖像做T=[9，9，2]的剛體變換。每一組測試中參考和浮動圖像按照文獻[14]的方法隨機添加20次強度不同的的萊斯噪聲。第四組測試的噪聲標(biāo)準(zhǔn)差為14，第五組的噪聲標(biāo)準(zhǔn)差為17，第六組的噪聲標(biāo)準(zhǔn)差為19。圖5為每組測試中20次配準(zhǔn)的成功率。

圖5 實驗二的配準(zhǔn)結(jié)果

由圖5可知，當(dāng)空間變換大并有噪聲存在時，NMI不能完成配準(zhǔn)任務(wù)，配準(zhǔn)率為0。第四組測試中本研究方法配準(zhǔn)成功率為100%，SCV為89.47%，CCRE為78.95%。第五組測試中本研究方法配準(zhǔn)成功率為89.45%，SCV為84.21%，CCRE為47.37%。第六組測試中本研究方法配準(zhǔn)成功率為89.47%，SCV為68.42%，CCRE為52.63%。由圖5可以看出，當(dāng)噪聲增大時，CCRE、SCV方法配準(zhǔn)成功率都會降低，本研究提出的方法仍可達到較高的配準(zhǔn)成功率。

表1列出以上六組測試中，成功配準(zhǔn)得到的空間參數(shù)誤差的平均值。

可以看出，本研究方法在保證配準(zhǔn)成功率的同時，具有和其他方法一致甚至明顯優(yōu)于其他方法的配準(zhǔn)精度。

表1配準(zhǔn)精度

Table 1 Registering accuracy

5 結(jié)束語

在多模醫(yī)學(xué)圖像配準(zhǔn)中，傳統(tǒng)相似性測度曲線存在局部極值，為空間變換參數(shù)尋優(yōu)帶來困難。針對此問題，本研究提出一種改進的SCV作為相似性測度。實驗表明，該方法與NMI、CCRE和原SCV相比有以下優(yōu)點：(1)當(dāng)參考圖像和浮動圖像的空間差別不大時，這四種相似性測度都能完成配準(zhǔn)任務(wù)，而本研究方法精度最高。(2)當(dāng)參考圖像和浮動圖像的空間差別較大或受到較大噪聲干擾時，本研究方法的配準(zhǔn)成功率最高。如前所述，本研究方法的相似性測度曲線較為光滑，從而提高配準(zhǔn)的魯棒性。(3)綜合比較，本研究方法在保證了配準(zhǔn)精度的同時對較大的空間變換和較強的噪聲具有魯棒性。該方法可以進一步擴展至3D醫(yī)學(xué)圖像配準(zhǔn)中。今后將進一步研究該方法的3D配準(zhǔn)應(yīng)用。