DBS術(shù)前圖像與ICBM-152圖譜的配準(zhǔn)算法

倪楊陽(yáng)，鄭慧芬，曹勝武，羅守華

1.東南大學(xué) 生物科學(xué)與醫(yī)學(xué)工程學(xué)院，江蘇 南京 210096；2.南京醫(yī)科大學(xué)附屬腦科醫(yī)院 老年神經(jīng)科，江蘇 南京 210029；3.江蘇省人民醫(yī)院 神經(jīng)外科，江蘇 南京 210029

引言

腦深部電刺激（Deep Brain Stimulation，DBS）手術(shù)，對(duì)于帕金森病等疾病是一種有效的治療手段。在帕金森病治療中，腦深部電刺激主要通過在基底核團(tuán)植入電極進(jìn)行刺激，從而達(dá)到治療效果?；颊咴诮邮蹹BS手術(shù)治療前，首先需要進(jìn)行手術(shù)計(jì)劃及術(shù)前導(dǎo)航。在局部麻醉?xiàng)l件下對(duì)患者安裝手術(shù)頭架，并進(jìn)行頭部MR以及CT掃描，通過影像融合確定立體定向頭架和電極植入點(diǎn)的定位。術(shù)前導(dǎo)航過程是一個(gè)直接決定術(shù)后效果的過程，對(duì)腦深部各個(gè)組織的識(shí)別至關(guān)重要，而人腦圖譜是反應(yīng)人類大腦解剖位置信息以及形狀信息的有力工具，兩者之間進(jìn)行配準(zhǔn)是一個(gè)值得研究的問題。

在DBS手術(shù)領(lǐng)域已經(jīng)有許多基于人腦圖譜的研究。Andreas等[1]將術(shù)前術(shù)后的影像數(shù)據(jù)同時(shí)向人腦圖譜配準(zhǔn)，通過圖譜中的特定標(biāo)記獲得丘腦底核的位置，并在此基礎(chǔ)上實(shí)現(xiàn)了DBS術(shù)后評(píng)估的工具包。Videen等[2]根據(jù)術(shù)后CT圖像中的電極金屬偽影定位電極中心，并將手術(shù)前后圖像與Mai圖譜配準(zhǔn)。這兩種方法在一定程度上可以解決丘腦底核分割以及電極定位問題，但是由于腦圖譜和患者腦部存在差異，因此核團(tuán)分割結(jié)果誤差難以預(yù)測(cè)。Silva等[3]將患者數(shù)據(jù)使用非剛性配準(zhǔn)方法配準(zhǔn)至標(biāo)準(zhǔn)圖譜，解決了患者與圖譜結(jié)構(gòu)上的差異帶來的誤差，但是在配準(zhǔn)過程中會(huì)導(dǎo)致電極的扭曲，會(huì)在結(jié)果中引入誤差。Castro等[4]比較了專家手動(dòng)定位與基于圖譜的算法以及基于非剛性配準(zhǔn)的算法并進(jìn)行了交叉驗(yàn)證。以上幾種方法都是將患者圖像配準(zhǔn)至圖譜，這樣做更方便使不同數(shù)據(jù)之間的結(jié)果量化，但是這樣做會(huì)導(dǎo)致病人的圖像由于非剛性配準(zhǔn)算法導(dǎo)致變形，從而使得算法結(jié)果無法應(yīng)用至術(shù)前的核團(tuán)識(shí)別中。

本文旨在通過使用DBS術(shù)前的臨床數(shù)據(jù)，對(duì)其進(jìn)行空間標(biāo)準(zhǔn)化后，將圖譜經(jīng)過粗配準(zhǔn)以及精確配準(zhǔn)兩個(gè)過程，實(shí)現(xiàn)臨床數(shù)據(jù)與腦圖譜數(shù)據(jù)的精確匹配。

1 算法設(shè)計(jì)

本文根據(jù)上述DBS手術(shù)影像數(shù)據(jù)的特點(diǎn)，選擇了目前使用較為廣泛的且有較多子圖譜應(yīng)用的ICBM-152圖譜作為標(biāo)準(zhǔn)圖譜數(shù)據(jù)。本文首先以患者圖像中的中矢狀面信息為依據(jù)，將患者數(shù)據(jù)標(biāo)準(zhǔn)化至MNI坐標(biāo)空間，隨后將標(biāo)準(zhǔn)的ICBM-152圖譜與標(biāo)準(zhǔn)空間下的T1數(shù)據(jù)經(jīng)過分段線性配準(zhǔn)以及非剛性配準(zhǔn)，得到病人的個(gè)性化圖譜，見圖1。

圖1 本文算法流程圖

1.1 患者圖像與腦圖譜的粗配準(zhǔn)

為實(shí)現(xiàn)臨床數(shù)據(jù)與圖譜的精確配準(zhǔn)，首先需要將兩者進(jìn)行粗略配準(zhǔn)，使二者在腦組織的大小，形狀以及外輪廓上趨于一致。為了完成這一目標(biāo)，首先需要將患者的腦圖像使用的坐標(biāo)空間轉(zhuǎn)換成與圖譜中人腦方向一致的MNI坐標(biāo)空間[5]。MNI坐標(biāo)空間的定義，見圖2，該坐標(biāo)空間以大腦前連合、后連合以及腦中矢狀面作為基準(zhǔn)，將前后連合線作為水平軸，將中矢狀面作為垂直面，該坐標(biāo)空間被神經(jīng)外科醫(yī)生廣泛認(rèn)可。

圖2 MNI坐標(biāo)空間示意圖

完成坐標(biāo)轉(zhuǎn)換首先需要識(shí)別中矢狀面，本文采用的是基于K-L距離提取中矢狀面的方法[6]。K-L距離衡量的是相同事件空間里的兩個(gè)概率分布的差異情況的一種度量。假設(shè)兩個(gè)隨機(jī)變量p={pi}和q={qi}，其中pi和qi表示第i種事件的發(fā)生概率。則K-L距離的定義為：

中矢狀面是經(jīng)過大腦半球間裂隙的平面，見圖3，中矢狀面上腦脊液灰度值出現(xiàn)的概率最大，因此與其他矢狀面的K-L距離最大。隨后將中矢狀面校正至垂直后需要將AP平面校正至水平面，同時(shí)將AC-PC旋轉(zhuǎn)至水平。

圖3 距離中矢狀面兩側(cè)2 cm的矢狀面與中矢狀面的比較圖

由于人腦與圖譜差異較大，因此需要對(duì)本文采用分段仿射變換將腦圖譜中的人腦大小變換成與患者較為接近的大小。借鑒Talairach-Nowinski標(biāo)識(shí)點(diǎn)[7]的定義，找出A、P、L、R、S、I這6個(gè)皮層標(biāo)識(shí)點(diǎn)[8]（圖4）。根據(jù)這些點(diǎn)計(jì)算出患者與圖譜的根據(jù)NMI空間定義的分塊的12個(gè)長(zhǎng)方體的每個(gè)維度的線性比例系數(shù)[9]。

圖4 AP、VAC、VPC上的標(biāo)識(shí)點(diǎn)

根據(jù)計(jì)算出的每個(gè)分塊對(duì)應(yīng)的每個(gè)維度上的縮放參數(shù)，對(duì)腦圖譜進(jìn)行分段縮放[10]，使圖像與圖譜在大小上達(dá)到一致。

1.2 圖譜與患者圖像間的非剛性配準(zhǔn)

在對(duì)圖譜完成上述變換后，圖譜與人腦的腦組織紋理仍然不匹配，因此需要進(jìn)行更為精細(xì)的非剛性配準(zhǔn)。本文所采用的是基于B樣條的自由形變模型的三維非剛性配準(zhǔn)[11]。

基于B樣條函數(shù)的配準(zhǔn)方法主要利用控制點(diǎn)和B樣條函數(shù)來描述非線性幾何變換。在三維空間中基于B樣條函數(shù)的自由形變可以描述為3個(gè)一維三次B樣條函數(shù)的張量積：

其中P表示控制點(diǎn)網(wǎng)格，βi表示i階B樣條函數(shù)，具體形式如下：

2 實(shí)驗(yàn)結(jié)果

本文實(shí)驗(yàn)算法設(shè)計(jì)及實(shí)驗(yàn)過程在Matlab以及Visual Studio 2010上完成，并使用了分割及配準(zhǔn)工具包[12]（Insight Segmentation and Registration Toolkit，ITK）以及視覺化工具函式庫(kù)[13]（Visualization ToolKi，VTK）中的部分算法。

對(duì)于B樣條自由形變配準(zhǔn)算法，ITK中使用了BSplineTransform類對(duì)其進(jìn)行了封裝，通過調(diào)用SetTransformDomainMeshSize函數(shù)，將每個(gè)維度使用的網(wǎng)格數(shù)設(shè)置為5個(gè)?？紤]網(wǎng)格點(diǎn)在3個(gè)維度上的移動(dòng)量，此時(shí)配準(zhǔn)所需要優(yōu)化的參數(shù)總數(shù)達(dá)到了375個(gè)。這個(gè)數(shù)量的優(yōu)化參數(shù)下，很難通過搜索類算法找到全局最優(yōu)值，因此本文采用了L-BFGS優(yōu)化算法[14]，該方法對(duì)于大量參數(shù)優(yōu)化有著很強(qiáng)的解決能力，其在ITK中被封裝為L(zhǎng)BFGSBOptimizerv 4類，本文中設(shè)置其代價(jià)函數(shù)收斂因子為1×107，迭代次數(shù)為200，漸變收斂容差為1×10-6。

2.1 原始數(shù)據(jù)

原始數(shù)據(jù)為7名帕金森患者，使用德國(guó)西門子公司3.0 T MRI系統(tǒng)，在快速自旋回波序列T1相重復(fù)時(shí)間（Time of Repetition，TR）為1600 ms，回波時(shí)間（Time of Echo，TE）為2.48 ms，可視范圍（Field of View，F(xiàn)OV）250 mm×250 mm，像素512×512，層厚1.5 mm條件下，或在快速自旋回波序列T2相TR 5000 ms，TE 84.00 ms，F(xiàn)OV 229 mm×229 mm，像素384×384，層厚2 mm條件下獲取的腦部MR數(shù)據(jù)。

術(shù)前數(shù)據(jù)為橫斷位的MR數(shù)據(jù)，圖譜數(shù)據(jù)為ICBM-152圖的標(biāo)準(zhǔn)數(shù)據(jù)（圖5）。

圖5 DBS術(shù)前臨床圖像以及ICBM-152圖譜部分橫斷位數(shù)據(jù)

2.2 圖譜與患者圖像配準(zhǔn)結(jié)果

分段線性配準(zhǔn)算法的結(jié)果，見圖6。通過網(wǎng)格細(xì)分比較，可以看出分段線性配準(zhǔn)前圖譜與圖像腦部差異較大，而進(jìn)行分段線性配準(zhǔn)后圖像與圖譜在結(jié)構(gòu)大小上趨于一致。

對(duì)分段線性變換后的圖譜進(jìn)行非剛性變換后的結(jié)果，見圖7。從圖中可以看出，圖譜的外輪廓不再是對(duì)稱的。同時(shí)，圖譜與圖像的內(nèi)部紋理更加趨于一致。

分別計(jì)算各個(gè)階段圖譜與臨床圖像的互信息值[15-16]，見表1。結(jié)合表1可以在分段仿射變換后互信息的提升了0.511，而在整個(gè)算法結(jié)束后互信息值提升了0.659。

圖6 分段線性配準(zhǔn)前后圖譜與患者M(jìn)R比較

圖7 非剛性配準(zhǔn)前后圖譜與患者M(jìn)R比較

表1 配準(zhǔn)各階段互信息值

3 總結(jié)

本文基于ICBM-152圖譜以及DBS術(shù)前圖像，首先通過中矢狀面的矯正，使圖像與圖譜空間達(dá)到一致，隨后使用分段線性配準(zhǔn)使圖像與圖譜的大腦大小達(dá)到一致，最后使用非剛性配準(zhǔn)算法使圖像與圖譜達(dá)到較為準(zhǔn)確的匹配。整個(gè)配準(zhǔn)流程都使用互信息作為相似性測(cè)度，適用于MR圖像不同的加權(quán)結(jié)果。通過計(jì)算各個(gè)配準(zhǔn)階段的二者的互信息值，發(fā)現(xiàn)最后的互信息值與初始的互信息值相比有一定程度的增加，配準(zhǔn)后的圖譜在輪廓以及內(nèi)部紋理上與臨床圖像達(dá)到了一定程度上的吻合。

本文使用的互信息測(cè)度具有較好的魯棒性，但運(yùn)算量較大，由于本文所采用的數(shù)據(jù)為三維MR序列，因此該算法處理所需的運(yùn)算時(shí)間較長(zhǎng)。其次，通過對(duì)比可以發(fā)現(xiàn)，配準(zhǔn)后的圖譜在局部仍有失配的現(xiàn)象。這兩個(gè)問題的解決是本文的方法用于實(shí)際應(yīng)用中的關(guān)鍵。