李碧草 李潤川 劉洲峰 李春雷 王宗敏

1(中原工學(xué)院電子信息學(xué)院 河南 鄭州 450007) 2(鄭州大學(xué)信息工程學(xué)院 河南 鄭州 450001) 3(鄭州大學(xué)互聯(lián)網(wǎng)醫(yī)療與健康服務(wù)河南省協(xié)同創(chuàng)新中心 河南 鄭州 450052)

0 引 言

在醫(yī)學(xué)成像、計(jì)算機(jī)視覺、遙感成像領(lǐng)域，圖像配準(zhǔn)是一個(gè)非常重要的技術(shù)。配準(zhǔn)算法在醫(yī)學(xué)圖像處理研究中，如計(jì)算機(jī)輔助診斷和腫瘤放射治療方面有著重要的作用。醫(yī)學(xué)成像過程中，由于病人體位的變化及器官和組織的移動(dòng)，不同時(shí)間采集的醫(yī)學(xué)圖像間往往存在空間變換，醫(yī)生為了更好地掌握病灶的信息，往往需要從不同模態(tài)的圖像(比如MRI、CT和US等醫(yī)學(xué)圖像)獲取病灶的互補(bǔ)信息，從而更精確地診斷病情。為了給醫(yī)生提供不同模態(tài)圖像的互補(bǔ)信息，需要一種魯棒的多模態(tài)配準(zhǔn)算法對圖像進(jìn)行配準(zhǔn)融合。

目前的圖像配準(zhǔn)方法通?？梢苑譃榛谔卣骱突诿芏葍纱箢怺1]?；谔卣鞯呐錅?zhǔn)算法包括特征提取和特征匹配，該類方法的配準(zhǔn)精度很大程度上取決于特征點(diǎn)的提取。由于人體組織和器官及醫(yī)學(xué)成像設(shè)備自身的特點(diǎn)，利用常見的特征提取方法提取特征點(diǎn)相對較難而且精度不高。本文研究基于密度的配準(zhǔn)方法。在基于密度的配準(zhǔn)方法中，基于互信息的配準(zhǔn)算法[2-3]自問世之后就備受關(guān)注。一些學(xué)者對基于互信息的配準(zhǔn)算法進(jìn)行改進(jìn)，提出歸一化的互信息配準(zhǔn)算法[4]、結(jié)合梯度信息的互信息配準(zhǔn)算法[5]、基于互累積剩余熵相似度的配準(zhǔn)算法[6]等。這些算法中的相似性測度都是由香農(nóng)熵構(gòu)造的，然而Antolin等[7]指出香農(nóng)熵的擴(kuò)展性使得基于香農(nóng)熵的信息論相似性測度并未考慮兩個(gè)獨(dú)立隨機(jī)變量間的相互作用。為此，學(xué)者引入一個(gè)具有非擴(kuò)展性的Tsallis熵，利用其性質(zhì)[8]提出一種相似性測度，以此相似度作為配準(zhǔn)標(biāo)準(zhǔn)進(jìn)行醫(yī)學(xué)圖像配準(zhǔn)。受文獻(xiàn)[7]的啟發(fā)，本文利用Arimoto非擴(kuò)展熵提出一種新的相似性測度，以此建立多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)模型，提高圖像配準(zhǔn)的精度。

1 非擴(kuò)展熵相似性測度

與香農(nóng)熵不同，Arimoto熵是一種非擴(kuò)展熵，也是香農(nóng)熵的一種推廣。本文利用Arimoto熵的性質(zhì)，構(gòu)造一種相似性測度。

1.1 Arimoto熵

假設(shè)一個(gè)離散隨機(jī)變量X，其概率分布為P=(p1,p2,…,pM)，則X的Arimoto熵定義為[9]：

(1)

依據(jù)洛必達(dá)法則求得當(dāng)α→1時(shí)，Arimoto熵的極限就等于香農(nóng)熵，因此Arimoto熵是一種廣義熵。Boekee等[10]研究了Arimoto熵的一些重要的性質(zhì)，下面僅討論本文用到的其中幾個(gè)性質(zhì)：

(1) 非負(fù)性。

Aα(X)≥0α>0,α≠1

(2)

(2) 非擴(kuò)展性。另給一個(gè)離散隨機(jī)變量Y，當(dāng)X與Y相互獨(dú)立時(shí)，它們的聯(lián)合Arimoto熵與邊緣Arimoto熵存在如下關(guān)系：

(3)

可見，除了兩個(gè)變量邊緣Arimoto熵的和之外，聯(lián)合Arimoto熵還包含一個(gè)乘積項(xiàng)，該乘積項(xiàng)表示兩個(gè)隨機(jī)變量間的相互作用。

(3) 凸性。

Aα(tX1+(1-t)X2)≥tAα(X1)+(1-t)Aα(X2)

(4)

(4) 對稱性。對稱性是指互換概率分布中任意兩個(gè)概率分量的位置，Arimoto熵的值不會(huì)發(fā)生變化，即：

Aα(…,pi,…,pj,…)=Aα(…,pj,…,pi,…)

(5)

此性質(zhì)與香農(nóng)熵的對稱性一樣，根據(jù)Arimoto熵的定義很容易驗(yàn)證。

(5) 極值性。

(6)

式(6)中等號(hào)成立的條件是：當(dāng)且僅當(dāng)pi=1/M，i=1,2,…,M。極值性也是一個(gè)很重要的性質(zhì)，根據(jù)極值性和非負(fù)性可以看到，當(dāng)M和α的值給定時(shí)，Arimoto熵的值是有界的。

1.2 Jensen-Arimoto散度

受文獻(xiàn)[11]的啟發(fā)，根據(jù)詹森香農(nóng)散度的定義并結(jié)合Arimoto熵的性質(zhì)，提出了一個(gè)新的散度測量——詹森Arimoto散度(Jensen-Arimoto Divergence，JAD)。與詹森香農(nóng)散度一樣，JAD也可以用來量化隨機(jī)變量之間的相似程度。

定義1假定X=(x1,x2,…,xM)是一個(gè)隨機(jī)變量，其概率分布為P=(p1,p2,…,pM)，則概率分布P的JAD定義為[12]：

α>0,α≠1

(7)

式中：Aα(·)代表Arimoto熵；ωi表示加權(quán)因子，而且ωi≥0，∑ωi=1。由于當(dāng)α趨向于1時(shí)Arimoto熵的極限是香農(nóng)熵，因此當(dāng)α→1時(shí)JAD的極限為詹森香農(nóng)散度(Jensen-Shannon Divergence，JSD)。

定理1當(dāng)p1,p2,…,pM為退化分布時(shí)，JAD能獲得最大值A(chǔ)α(ω)，這里pi=δij，δij表示克羅內(nèi)克符號(hào)(Kronecker Symbol)，即當(dāng)i=j時(shí)，δij=1，否則δij=0。

2 圖像配準(zhǔn)算法

在醫(yī)學(xué)成像領(lǐng)域，模態(tài)是指不同類型的成像方式，如磁共振成像、CT成像、超聲成像等。多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)中，待配準(zhǔn)的對象是不同模態(tài)的圖像。首先建立配準(zhǔn)模型，然后選擇合適的空間變換，再利用上文介紹的相似度作為目標(biāo)函數(shù)，最后采用擬牛頓法對目標(biāo)函數(shù)進(jìn)行優(yōu)化，獲得最終的變換參數(shù)，完成多模態(tài)醫(yī)學(xué)圖像的精確配準(zhǔn)。

2.1 圖像配準(zhǔn)框架

給定兩幅待配準(zhǔn)圖像，一幅記為參考圖像R，另一幅為浮動(dòng)圖像F。圖像配準(zhǔn)可以看作尋找待配準(zhǔn)圖像間的空間變換的過程。設(shè)x=(x,y,z)T為參考圖像上的任意一個(gè)像素點(diǎn)，它經(jīng)過空間變換后對應(yīng)于浮動(dòng)圖像中的點(diǎn)x′=(x′,y′,z′)T，兩點(diǎn)之間的關(guān)系表述為：

x′=Tμ(x)

(8)

式中:Tμ(·)是空間變換函數(shù)；μ為變換參數(shù)，這里的變換可以是剛體變換、仿射變換甚至是彈性形變。圖1表示兩幅圖像對應(yīng)坐標(biāo)之間的空間變換關(guān)系,(a)中圓點(diǎn)的坐標(biāo)為x，(b)中星號(hào)點(diǎn)的坐標(biāo)為x′。

(a) (b)圖1 圖像坐標(biāo)的空間變換關(guān)系

由此F與R的配準(zhǔn)過程可以看成如下優(yōu)化問題：

(9)

式中：S用來衡量兩幅圖像的相似性測度，當(dāng)R(x)和F(Tμ(x))完全配準(zhǔn)時(shí)，該相似度取得最大值，符號(hào)“°”代表空間變換Tμ作用于浮動(dòng)圖像。圖像配準(zhǔn)框架的示意圖如圖2所示。

圖2 圖像配準(zhǔn)框架示意圖

2.2 空間變換模型

由于本文的實(shí)驗(yàn)對象是腦部圖像，因此我們主要考慮剛體變換。剛體變換包含平移和旋轉(zhuǎn)，對于三維圖像來說，剛體變換有6個(gè)自由度：x、y、z方向上的平移量與繞x、y、z軸的旋轉(zhuǎn)角度。三維圖像的剛體變換可用矩陣形式表示為：

X′=RX+t

(10)

式中：X=(x,y,z)T和X′=(x′,y′,z′)T分別代表變換前后圖像中像素點(diǎn)的坐標(biāo)；R是3×3的旋轉(zhuǎn)矩陣；t為3×1的平移向量。旋轉(zhuǎn)矩陣R可以用三個(gè)矩陣的乘積來表示：

R=r(x)r(y)r(z)

(11)

矩陣r(x)、r(y)和r(z)可以用三個(gè)旋轉(zhuǎn)角度表示：

(12)

(13)

(14)

式中：α、β、γ分別代表圖像繞x、y、z軸的旋轉(zhuǎn)角度。由此三維剛體圖像配準(zhǔn)的目的就是通過某種優(yōu)化方法尋找使得待配準(zhǔn)圖像間的相似度最大時(shí)的6個(gè)變換參數(shù)。

2.3 相似性測度

本文利用第1節(jié)中構(gòu)造的基于Arimoto熵的相似度來衡量兩幅圖像間的相似程度，將兩幅圖像看作隨機(jī)變量，那么用它們的灰度概率分布計(jì)算其相似度。給定參考圖像R和浮動(dòng)圖像F，及它們之間的空間變換Tμ，假設(shè)用于估計(jì)概率分布的箱子數(shù)(Number of bins)為M，讓f=(f1,f2,…,fM)和r=(r1,r2,…,rM)分別為F(Tμ(x))和R(x)的灰度級(jí)。式(7)中，令JAD定義中的pi=pi(F(Tμ(x))|R(x))，i=1,2,…,M，表示已知參考圖像前提下變換的浮動(dòng)圖像的條件概率分布，使用兩幅圖像的灰度級(jí)可以將條件概率分布重新寫為pi=p(F=fj|R=ri)=p(fj|ri)，然后再用參考圖像的概率分布作為加權(quán)因子，即ωi=p(R(x))=p(R=ri)=p(ri)，將pi和ωi代入JAD可得：

JAα(F(Tμ(x)),R(x))=

(15)

將式(15)中的JA測度替換S代入式(9)的配準(zhǔn)框架，接下來的工作就變成如何求該配準(zhǔn)框架最優(yōu)解的問題。

2.4 優(yōu)化方案

與梯度下降的優(yōu)化方法相比，牛頓法中二階信息的應(yīng)用能夠提供更好的理論收斂性質(zhì)[13]，但是海森矩陣(Hessianmatrix)及其求逆的計(jì)算非常耗時(shí)，為此有學(xué)者提出擬牛頓法。該方法不需要計(jì)算海森矩陣的逆，而是用它的估計(jì)來代替，因此并不需要計(jì)算目標(biāo)函數(shù)的二階導(dǎo)數(shù)。本文采用擬牛頓法的一個(gè)變種——L-BFGS優(yōu)化方法[14]，來尋找目標(biāo)函數(shù)的最大值。利用泰勒級(jí)數(shù)對目標(biāo)函數(shù)進(jìn)行展開，只保留前三項(xiàng)得到：

S(μ+Δμ)≈S(μ)+ΔμT·▽S(μ)+

(16)

式中：Δμ表示參數(shù)向量的增量。對于L-BFGS方法，參數(shù)向量第k+1步迭代更新可以由第k步得到：

μ(k+1)=μ(k)-(H(k))-1·▽S(μ(k))

(17)

式中：(H(k))-1代表每次迭代的步長，是海森矩陣逆的估計(jì)，并非真實(shí)計(jì)算的海森矩陣的逆；▽表示梯度算子。接下來需要計(jì)算目標(biāo)函數(shù)關(guān)于變換參數(shù)的導(dǎo)數(shù)：

(18)

式中：n表示空間變換模型中參數(shù)(自由度)的個(gè)數(shù)。對于三維圖像的仿射變換，n=12。

要計(jì)算目標(biāo)函數(shù)關(guān)于參數(shù)μ的導(dǎo)數(shù)，首先要推導(dǎo)出連續(xù)的目標(biāo)函數(shù)表達(dá)式，本文將采用基于B樣條的Parzen窗來估計(jì)圖像的概率密度函數(shù)。設(shè)β(0)和β(3)分別代表零階和三階B樣條函數(shù)，由于參考圖像在配準(zhǔn)過程中并不受變換參數(shù)μ的影響，因此采用零階B樣條函數(shù)來估計(jì)參考圖像的概率分布。對于空間變換后的浮動(dòng)圖像，三階B樣條函數(shù)用于估計(jì)其概率密度函數(shù)。那么根據(jù)多維B樣條函數(shù)的可分離性，R(x)和F(Tμ(x))的聯(lián)合概率密度函數(shù)可由它們各自的邊緣概率密度函數(shù)相乘得到：

(19)

通過計(jì)算得到相似性測度S關(guān)于變換參數(shù)μ的導(dǎo)數(shù)如下：

(20)

(21)

式中:β′(3)表示三階B樣條函數(shù)的導(dǎo)數(shù)；?F(t)/?t代表變換的浮動(dòng)圖像F(Tμ(x))的梯度；?(Tμ(x))/?μ為空間變換模型關(guān)于參數(shù)μ的導(dǎo)數(shù)。

本文的配準(zhǔn)實(shí)驗(yàn)中，采用的L-BFGS優(yōu)化算法的停止條件為：相鄰兩步迭代的目標(biāo)函數(shù)值之差小于0.01或者迭代次數(shù)大于預(yù)先設(shè)定的值，本文設(shè)定該值為100。

3 實(shí)驗(yàn)結(jié)果與分析

為了驗(yàn)證本文算法的性能，我們選擇不同模態(tài)的三維頭部磁共振圖像進(jìn)行實(shí)驗(yàn)。這些測試圖像來源于蒙特利爾神經(jīng)醫(yī)學(xué)部的BrainWeb數(shù)據(jù)庫[15]，包含三種不同類別的磁共振圖像：T1加權(quán)、T2加權(quán)和PD加權(quán)的MRI。三組數(shù)據(jù)的尺寸都為181×217×181，在物理坐標(biāo)下每個(gè)體素的大小為1 mm×1 mm×1 mm。而且這三組圖像的所有層都是相互對應(yīng)的。圖3顯示了三個(gè)MRI圖像的三個(gè)正交層。

(a) T1加權(quán)MR

(b) T2加權(quán)MR

(c) PD加權(quán)MR，每一行分別表示軸向面、矢狀面和冠狀面圖3 腦部磁共振圖像

接下來的實(shí)驗(yàn)采用圖3中的(a)、(b)和(c)作為實(shí)驗(yàn)圖像來執(zhí)行配準(zhǔn)算法。使用本文算法進(jìn)行圖像配準(zhǔn)時(shí)，需要考慮兩個(gè)主要因素：參數(shù)α對Arimoto熵的影響和計(jì)算熵時(shí)所選取的圖像塊尺寸。

為了驗(yàn)證所提算法的性能，三種MR圖像組成三對測試圖像：MR T1&MR T2、MR T1&MR PD及MR PD&MR T2，將每組測試數(shù)據(jù)中的前者作為參考圖像，然后利用已知的3D剛體變換(真實(shí)值)對后者執(zhí)行變換，把變換后的圖像作為浮動(dòng)圖像，如圖4所示。

(a) MR T1圖像作為參考圖像 (b) MR T2圖像 (c) 3D剛體變換后的MR T2圖像作為浮動(dòng)圖像圖4 測試圖像

在每組測試圖像中，選擇30個(gè)3D剛體變換對后者圖像進(jìn)行變換生成30幅浮動(dòng)圖像，這些變換的平移及旋轉(zhuǎn)參數(shù)隨機(jī)產(chǎn)生，并服從以(10 mm,8°)為均值，(2 mm,2°)為標(biāo)準(zhǔn)差的正態(tài)分布。由于三組測試圖像最初是對齊的，所以這30個(gè)剛體變換的參數(shù)即為配準(zhǔn)的真實(shí)值。為了能定量地評(píng)估配準(zhǔn)方法的精度，我們將利用配準(zhǔn)方法得到的參數(shù)值與真實(shí)值的差異作為配準(zhǔn)誤差，然后比較所提算法(JAD)與另外兩種方法NMI和CCRE的配準(zhǔn)誤差。

表1展示了這三種方法應(yīng)用到三組實(shí)驗(yàn)數(shù)據(jù)的配準(zhǔn)誤差?？梢钥闯?，在處理多模態(tài)圖像配準(zhǔn)時(shí)，與另外兩種方法相比本文方法獲得了較低的配準(zhǔn)誤差。三組實(shí)驗(yàn)的配準(zhǔn)精度都達(dá)到了亞像素級(jí)，而基于CCRE相似度的算法配準(zhǔn)效果相對較差，特別在配準(zhǔn)T2加權(quán)MR與PD加權(quán)MR時(shí)，誤差超過了一個(gè)像素。收斂范圍也是衡量一個(gè)相似度優(yōu)劣的重要指標(biāo)，也就是說，當(dāng)待配準(zhǔn)圖像之間存在不同大小的變換時(shí)，尤其是當(dāng)圖像間的誤匹配很大時(shí)，一個(gè)魯棒的配準(zhǔn)算法也能夠獲得最優(yōu)變換。為了測試JAD相似度的收斂范圍，使用三維的MR T1和MR PD作為測試圖像，并執(zhí)行五組實(shí)驗(yàn)，每一組中采用30個(gè)3D剛體變換生成30幅浮動(dòng)圖像。

表1 執(zhí)行30次3D剛體配準(zhǔn)實(shí)驗(yàn)的誤差均值±標(biāo)準(zhǔn)差

五組實(shí)驗(yàn)中采用的剛體變換參數(shù)分別來自五個(gè)不同大小的區(qū)間：[-10,10]，[-20,20]，[-30,30]，[-40,40]和[-50,50]，單位是mm或度數(shù)。為定量分析配準(zhǔn)結(jié)果，引入成功率的概念，算法的成功率定義為在所有的實(shí)驗(yàn)中成功次數(shù)所占的比例，所謂的“成功”代表優(yōu)化算法在100次迭代中能夠收斂于最優(yōu)解。圖5描述了利用三種方法執(zhí)行五組配準(zhǔn)實(shí)驗(yàn)的成功率?？梢钥闯?，當(dāng)初始變換較小時(shí)，三種方法都能獲得到很高的成功率，而隨著變換的增大，CCRE方法的配準(zhǔn)成功率下降很明顯，JAD和NMI相似度的成功率也有所降低，但JAD能夠提供比NMI稍高的成功率。綜上所述，JAD擁有較大的收斂范圍，并且當(dāng)圖像間存在較大變換時(shí)，可以得到更多的配準(zhǔn)次數(shù)。

圖5 不同大小變換下三種配準(zhǔn)方法的成功率比較

4 結(jié) 語

針對傳統(tǒng)基于互信息相似度的圖像配準(zhǔn)方法，本文提出一種基于非擴(kuò)展熵相似度的多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)算法。首先，研究Arimoto非擴(kuò)展熵的性質(zhì)，利用Arimoto熵構(gòu)造相似性測度，用其衡量兩幅圖像間的相似度；其次，結(jié)合空間變換模型建立圖像配準(zhǔn)框架；然后，利用基于B樣條的帕曾窗方法估計(jì)待配準(zhǔn)圖像間的聯(lián)合概率分布，得到連續(xù)的目標(biāo)函數(shù)；最后，采用擬牛頓優(yōu)化方法對配準(zhǔn)模型進(jìn)行求解，得到待配準(zhǔn)圖像間最優(yōu)的空間變換參數(shù)，實(shí)現(xiàn)多模態(tài)醫(yī)學(xué)圖像的精確配準(zhǔn)。三維腦部磁共振圖像配準(zhǔn)實(shí)驗(yàn)結(jié)果表明，與歸一化互信息算法和基于互累積剩余熵的配準(zhǔn)方法相比，本文算法能夠獲得更高的配準(zhǔn)精度。未來研究將考慮三維彈性配準(zhǔn)，并進(jìn)一步提高配準(zhǔn)精度。