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

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

1 引 言

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

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

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

2 理論與方法2.1 SCV原理

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

(1)

(2)

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

2.2 改進(jìn)的插值方法

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

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

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

(3)

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

wi=f(dix)·f(diy)

(4)

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

圖1 改進(jìn)的插值方法2D示意圖

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

(1)f(x)≥0；

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

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

(5)

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

(6)

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

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

(7)

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

2.3 本研究的相似性測度

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

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

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

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

(8)

(9)

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

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

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

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

4 多模圖像配準(zhǔn)實(shí)驗(yàn)

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

(10)

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

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

由上圖可知，第一組測試中，當(dāng)空間變換不大時(shí)， NMI、CCRE、SCV和本研究的方法成功率均達(dá)到了100%。第二組測試中，當(dāng)空間變換加大時(shí)，本研究方法配準(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)浮動(dòng)圖像相比參考圖像的空間變化增加時(shí)，NMI、CCRE、SCV方法配準(zhǔn)成功率都會(huì)降低，本研究提出的方法仍然可以達(dá)到較高的配準(zhǔn)成功率。

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

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

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

由圖5可知，當(dāng)空間變換大并有噪聲存在時(shí)，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)噪聲增大時(shí)，CCRE、SCV方法配準(zhǔn)成功率都會(huì)降低，本研究提出的方法仍可達(dá)到較高的配準(zhǔn)成功率。

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

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

表1配準(zhǔn)精度

Table 1 Registering accuracy

5 結(jié)束語

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