DBS術前圖像與ICBM-152圖譜的配準算法

倪楊陽，鄭慧芬，曹勝武，羅守華

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

引言

腦深部電刺激（Deep Brain Stimulation，DBS）手術，對于帕金森病等疾病是一種有效的治療手段。在帕金森病治療中，腦深部電刺激主要通過在基底核團植入電極進行刺激，從而達到治療效果。患者在接受DBS手術治療前，首先需要進行手術計劃及術前導航。在局部麻醉條件下對患者安裝手術頭架，并進行頭部MR以及CT掃描，通過影像融合確定立體定向頭架和電極植入點的定位。術前導航過程是一個直接決定術后效果的過程，對腦深部各個組織的識別至關重要，而人腦圖譜是反應人類大腦解剖位置信息以及形狀信息的有力工具，兩者之間進行配準是一個值得研究的問題。

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

本文旨在通過使用DBS術前的臨床數據，對其進行空間標準化后，將圖譜經過粗配準以及精確配準兩個過程，實現臨床數據與腦圖譜數據的精確匹配。

1 算法設計

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

圖1 本文算法流程圖

1.1 患者圖像與腦圖譜的粗配準

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

圖2 MNI坐標空間示意圖

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

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

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

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

圖4 AP、VAC、VPC上的標識點

根據計算出的每個分塊對應的每個維度上的縮放參數，對腦圖譜進行分段縮放[10]，使圖像與圖譜在大小上達到一致。

1.2 圖譜與患者圖像間的非剛性配準

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

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

其中P表示控制點網格，βi表示i階B樣條函數，具體形式如下：

2 實驗結果

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

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

2.1 原始數據

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

術前數據為橫斷位的MR數據，圖譜數據為ICBM-152圖的標準數據（圖5）。

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

2.2 圖譜與患者圖像配準結果

分段線性配準算法的結果，見圖6。通過網格細分比較，可以看出分段線性配準前圖譜與圖像腦部差異較大，而進行分段線性配準后圖像與圖譜在結構大小上趨于一致。

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

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

圖6 分段線性配準前后圖譜與患者MR比較

圖7 非剛性配準前后圖譜與患者MR比較

表1 配準各階段互信息值

3 總結

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

本文使用的互信息測度具有較好的魯棒性，但運算量較大，由于本文所采用的數據為三維MR序列，因此該算法處理所需的運算時間較長。其次，通過對比可以發現，配準后的圖譜在局部仍有失配的現象。這兩個問題的解決是本文的方法用于實際應用中的關鍵。