腦深部電刺激（DBS）手術是中、晚期帕金森病患者的重要治療手段，電極在相應核團位置的精確度直接決定著手術的治療效果。目前單一影像學手段無法得到具有電極、核團及其位置關系的圖像，且由于丘腦底核體積小，邊界不明顯，無法對其直接進行分割。本文利用 DBS 手術特有的 T1、T2 及 SWI 加權方式的磁共振（MR）數據構建了完整的端到端 DBS 術后效果評估管線。首先將術前、術后患者圖像配準，標準化至同一坐標空間；隨后利用腦圖譜與術前數據的非剛性配準得到個性化的患者圖譜，從而得到術前的核團預測位置；之后再次利用術后圖像中的電極路徑及結合核團分割結果，得到電極和核團位置關系的三維圖像，有助于 DBS 術后效果的評估，為術后程控提供有效的信息。經分析，本文算法可以實現患者 DBS 手術圖像與腦圖譜間的配準，算法與專家評估的丘腦底核重心物理坐標的誤差為（1.590 ± 1.063）mm。本文解決了 DBS 手術電極植入位置的術后評估問題。
引用本文： 羅守華, 倪楊陽, 鄭慧芬, 曹勝武. 多模影像的腦深部電刺激術后效果評估管線的設計與實現. 生物醫學工程學雜志, 2019, 36(3): 356-363. doi: 10.7507/1001-5515.201711055 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
腦深部電刺激（deep brain stimulation，DBS）手術，對于帕金森病、特發性震顫等疾病是一種有效的治療手段。在帕金森病治療中，DBS 主要通過在丘腦底核植入電極進行刺激，從而達到恢復患者運動功能的效果。患者在接受各項評估并達到手術要求后可以進行 DBS 手術治療。其第一步是進行手術計劃，患者將在局部麻醉條件下安裝手術頭架，并進行頭部磁共振掃描（magnetic resonance imaging，MRI）以及計算機斷層掃描（computed tomography，CT），通過影像融合確定立體定向頭架和電極植入點的定位。隨后是進行植入手術，根據手術計劃在患者頭皮標記點位置鉆圓孔，將探測電極按照計劃路徑植入腦內，術中通過刺激器測試電生理信號來進行微調。這一過程中，電極植入的精確度直接關系著手術的效果。最后是進行術后程控，在此過程中出現的失語、認知障礙以及疼痛等副作用也與電極植入精確度有關。因此，術后對于電極植入位置的準確評估是一個迫切需要解決的問題。
目前，臨床上一般使用高分辨 MRI 獲得患者術前的全腦信息以及核團位置。在手術后，由于電極為金屬材質，不能使用高分辨 MRI 對植入部位進行成像，因此一般使用 CT 或低分辨 MRI 對患者頭部進行掃描。在臨床中，通常采用 T2 加權方式對丘腦底核進行成像。由于丘腦底核體積較小，且與周圍組織沒有明顯的邊界，因此直接通過觀察判斷其位置會引入誤差。其次，術后的影像由于分辨率較低，只能在其中觀察到電極位置。同時，患者在接受術前術后掃描時無法保證體位的一致。這三方面的因素為準確判斷術后電極與核團的相對位置帶來了困難。在實際操作中，T2 加權成像方式需要的時間較長，往往只能得到核團區域的腦部圖像，這對后續判斷以及算法的設計提出了更高的要求。
DBS 術后效果評估的主要難點為患者丘腦底核的精確分割，人腦圖譜在核團分割過程中起到很大作用。Horn 等[1]將術前術后的影像數據同時向人腦圖譜配準，通過圖譜中的特定標記獲得丘腦底核的位置，并在此基礎上實現了 DBS 術后評估的工具包。張芷齊等[2]在標準腦圖譜空間下將雙側丘腦底核刺激治療患者的高分辨率術前圖像與低分辨率的術后圖像進行配準，并重建 DBS 電極的路徑，為醫生的術后程控提供輔助指導。Videen 等[3]根據術后 CT 圖像中的電極金屬偽影定位電極中心，并將手術前后圖像與 Mai 圖譜配準。上述方法在一定程度上可以解決丘腦底核分割以及電極定位問題，但是由于腦圖譜和患者腦部存在差異，因此核團分割結果的誤差難以預測。da Silva 等[4]將患者數據使用非剛性配準方法配準至標準圖譜，解決了患者與圖譜結構上的差異帶來的誤差，但是在配準過程中會導致電極的扭曲，會在結果中引入誤差。郝斌[5]對丘腦底核進行手動標記，實現了術后丘腦底核區域可視化，并從術后 MRI 數據中得到電極位置，分析電極植入位置與 DBS 手術療效之間的關系。本文使用非剛性配準方法，將腦圖譜轉化為與患者匹配度較高的個性化圖譜，減少了術后數據到圖譜的配準過程，從而減少在此過程中引入的誤差，同時根據術前術后的數據特點設計了數據處理流程，實現了一個端到端 DBS 術后效果評估管線，即一套將 DBS 手術圖像用于術后效果評估的標準流程。
本文旨在利用帕金森病 DBS 術前高分辨圖像中的核團信息及術后低分辨圖像中的電極信息，將二者進行配準，從而實現電極與核團空間位置的可視化，為術后程控提供可靠信息。同時，在研究中引入標準腦圖譜，將圖譜與術前圖像配準，解決術前影像中核團難以分割的問題。在實現術前影像中丘腦底核分割后，將術前與術后的圖像配準，從而得到電極與核團的相互關系。
1 算法設計
本文根據上述 DBS 手術影像數據的特點，在確保充分利用信息的情況下設計了如圖 1 所示的影像評估管線，首先將術前 T1 以及 T2 數據配準至同一空間，同時將術前 T1 數據以及配準后的 T2 數據標準化至 MNI 坐標空間，得到標準空間下的 T1 以及 T2 數據。隨后將標準的 ICBM152 圖譜與標準空間下的 T1 數據經過分段線性配準以及非剛性配準，得到患者的個性化圖譜。最后為了避免手術導致的大腦形變引入誤差，使用非剛性配準算法將患者術后 T1 數據與術前的個性化圖譜進行配準，根據圖譜中的標記確定術后核團的位置，并且將圖像中的電極提取出來，從而分析二者的相對位置關系。
圖1
影像評估管線流程圖
Figure1.
Flow chart of image evaluation pipeline
1.1 圖像預處理
首先，對患者的 MRI 圖像進行配準，將術前 T2 加權的 MRI 圖像以及術后 T1 加權的 MRI 圖像分別與術前 T1 加權的 MRI 圖像使用分割及配準工具包（Insight Segmentation and Registration Toolkit，ITK）[6]中的剛性配準算法進行配準，將手術前后的數據標準化到同一坐標空間中。在剛性配準過程中，由于組織的邊緣以及灰度在不同的加權方式下無法保持一致，因此采用互信息作為圖像的相似性度量。互信息是信息論中的一個基本概念，通常用于描述兩個系統間的統計相關性，可以用熵來表達[7]。
系統 A 的熵定義為
 |
兩個系統的聯合熵為
 |
兩個系統間的互信息為
 |
當兩幅圖像空間位置完全一致時，可以用其中一幅圖像表達另一幅圖像的信息，可以使用聯合概率分布以及完全獨立的概率分布間的廣義距離估計互信息：
 |
1.2 患者圖像與圖譜配準
在完成配準過程后，所有患者圖像都被標準化到同一坐標空間內。由于核團在患者 MR 影像中不明顯，需要通過將患者 T2 加權圖像與 MNI-ICBM152 圖譜配準實現對核團區域的識別。為了完成這一目標，首先需要將患者的腦圖像使用的坐標空間轉換成與圖譜中人腦方向一致的 MNI 坐標空間[8]。MNI 坐標空間以大腦前聯合（anterior commissure，AC）、后聯合（posterior commissure，PC）和腦中隔面為基準，與過 AC-PC 的連線的水平面（簡稱 AP 平面），以及過 AC 和 PC 且垂直于腦中隔和水平面的兩個平面（簡稱 VAC 平面及 VPC 平面）組成的平面將大腦分為 12 個長方體區域，如圖 2 所示。
圖2
MNI 坐標空間示意圖
Figure2.
MNI coordinate space diagram
完成坐標轉換需要識別腦中隔面，并將其旋轉至垂直位置，本文采用的是 Volkau 等[9]提出的使用 K-L 距離提取腦中隔面的方法。該方法是一種基于統計學的三維方法，魯棒性較好。K-L 距離衡量的是相同事件空間里的兩個概率分布的差異情況。設兩個隨機變量 p = {pi}和 q = {qi}，其中 pi 和 qi 表示第 i 種事件的發生概率。K-L 距離的定義為：
 |
腦中隔是通過大腦半球間裂的平面，是含有腦脊液最多的矢狀面。如圖 3 所示，中隔面上腦脊液灰度值出現的概率最大，K-L 距離反映了兩個切片上灰度值出現概率的變化情況，由于中隔面的變化最大，因此與其他矢狀面的 K-L 距離最大。
圖3
距離中隔面兩側 2 cm、1 cm 的矢狀面與中隔面的比較圖
Figure3.
Comparison of the mid-sagittal plane and the two sides of the septum 2 cm, 1 cm from the plane
將中隔面校正至垂直后需要將 AP 平面校正至水平面，根據中隔面上的 AC、PC 點計算出需要圍繞旋轉的單位向量（nx，ny，nz）以及角度 θ，根據羅德里格斯旋轉公式，旋轉矩陣為：
 |
由于人腦的個體差異較大，若直接與腦圖譜執行精確配準算法，容易導致算法在優化時陷入局部極值影響最終結果，因此需要采用分段仿射變換，將腦圖譜中的人腦大小變換成與患者較為接近的大小，為后續精確配準算法提供一個良好的初值。使用 Liu 等[10]提出的方法確定 AC、PC 點。根據 MNI 標準空間定義，如圖 4 所示，在 AP 平面上找出：過中隔面的額葉最前點 A；過中隔面的枕葉最后點 P；過 VPC 的左臨時葉最邊緣點 L；過 VPC 的右臨時葉最邊緣點 R。在 VAC 平面上找出過中隔面的腦組織最下點 I，在 VPC 上找出過中隔面的頂葉最上點 S。根據這些點計算出患者與圖譜根據 MNI 空間定義的 12 個分塊在每個維度上的線性比例系數[11]。根據所得系數對腦圖譜進行分段縮放，隨后計算縮放后圖譜與患者 MR 圖像中的 AC、PC 點位置的差值，對整個圖譜進行平移，使圖譜與患者前后聯合重合。
圖4
AP、VAC、VPC 上的標識點
Figure4.
The identification points on the AP, VAC and VPC
1.3 圖譜與患者圖像間的非剛性配準
在對圖譜完成上述變換后，模擬人腦與患者腦部的大小基本達到了一致，但還存在局部失配的現象，因此需要進行更為精細的非剛性配準。非剛性配準通過使標準圖譜形變擬合患者圖像，使差異較大的不同個體圖像間達到配準。本文所采用的是基于 B 樣條的自由形變模型（free form deform，FFD）的三維非剛性配準[12]。
基于 B 樣條函數的配準方法主要利用控制點和 B 樣條函數來描述非線性幾何變換。移動位于樣條函數上的控制點，使樣條函數彎曲，樣條函數之間的空間也會產生扭曲，達到自由形變的效果。在三維空間中基于 B 樣條函數的自由形變可以描述為三個一維三次 B 樣條函數的張量積：
 |
其中 P 表示控制點網格，βi 表示 i 階 B 樣條函數，具體形式如下
 |
由于圖譜與患者 MR 圖像間的配準屬于多模態圖像配準，因此同樣采用互信息作為相似性測度。
1.4 手術前后核團位置評估
標準的人腦圖譜根據以上變換轉化為與患者匹配度更高的個性化圖譜，ICBM-152 圖譜包含了如圖 5 所示的丘腦底核精確分割信息，因此將計算所得的形變網格應用至核團標識像中，便可獲得患者核團位置的評估結果。
圖5
丘腦底核在圖譜中的標記
Figure5.
Mask of the subthalamic nucleus in the atlas
術后患者影像雖然在圖像預處理時標準化到了術前患者圖像坐標空間中，但是由于電極植入過程需要打開患者顱腔，顱內壓力變化以及隨后的電極植入過程都會導致大腦形變[13]，因此需要對術前圖像向術后圖像進行非剛性配準。
2 實驗結果
2.1 實驗對象
受試者為 4 名帕金森患者，均已知曉并同意本研究使用其臨床數據。受試者在江蘇省人民醫院完成雙側丘腦底核 DBS 手術，植入電極為美國美敦力公司 3389 電極。手術前，使用德國西門子公司 3.0T MRI 系統，在快速自旋回波序列 T1 相重復時間（time of repetition，TR）1 600 ms、回波時間（time of echo，TE）2.48 ms、可視范圍（field of view，FOV）250 mm × 250 mm、像素 512 × 512 以及層厚 1.5 mm 條件下，以矢狀位無間隔掃描獲取全腦數據，隨后在快速自旋回波序列 T2 相 TR 5 000 ms、TE 84.00 ms、FOV 229 mm × 229 mm、像素 384 × 384、層厚 2 mm 條件下以橫斷位掃描獲取核團部位數據。手術后，使用美國通用電氣公司 1.5T MRI 系統，在快速自旋回波 T1 相 TR 4 000 ms、TE 105.62 ms、FOV 250 mm × 20 mm、像素 512 × 512、層厚 2 mm 條件下以橫斷位掃描獲得帶有電極的 MRI 圖像。
本文采用基于蒙特利爾神經科學研究所標準立體空間標準圖譜 MNI-ICBM152。該圖譜是將 152 名志愿者的高空間分辨率掃描結果，通過仿射變換與 MNI305 進行對應，再將 152 個大腦進行平均所得到的標準模板，已經被國際大腦映射聯盟（International Consortium for Brain Mapping，ICBM）作為標準模板。本文實驗算法設計及實驗過程在 Matlab 以及 Visual Studio 2010 上完成，并使用 ITK 中的部分算法。
2.2 原始數據
原始數據如圖 6 所示。術前數據為用于觀察大腦全局形態的矢狀位 T1 全腦數據和用于觀察核團的 T2 橫斷位數據（由于儀器條件所限，T2 加權方式掃描過程較長，因此只提供核團部位的數據）。為了防止金屬電極對患者造成傷害，術后數據只提供核團部位的橫斷位低分辨率數據。
圖6
原始臨床數據
Figure6.
The original clinical data
2.3 圖像預處理結果
剛性配準結果如圖 7 所示。可見算法執行后，術前術后核團部位與術前全腦數據完好地配準到了同一坐標空間中。
圖7
不同圖像與術前 T1 加權 MR 剛性配準結果
Figure7.
Rigid registration results of different images with preoperative T1-weighted MR
中隔面校正算法的執行結果如圖 8 所示，可見該算法較好地將中隔面校正到了與 x 軸垂直的狀態。
圖8
中隔面校正算法執行前（左側）后（右側）
Figure8.
Images before (left) and after (right) the middle septal correction
圖 9 為將圖像序列按照 AC-PC 點旋轉的前后對比效果，其中右側圖像中的患者大腦已校正至 MNI 標準空間。
圖9
AP 平面校正執行前（左側）后（右側）
Figure9.
Images before (left) and after (right) the AP plane correction
2.4 圖譜與患者圖像配準結果
由各組圖像的配準結果可以看出，ICBM152 圖譜的腦部較為狹長，將配準前后的圖譜與圖像進行融合，可得到如圖 10 所示的分塊融合結果。從中可以看到，經過配準的圖譜與臨床數據的腦組織在大小上更加一致。為了對分段仿射配準的效果進行量化，本文計算了配準前后圖譜與圖像之間的歸一化互信息值，其結果如表 1 所示。歸一化互信息值的區間為[1，2]，互信息值越大代表兩幅圖像間的配準效果越好。
圖10
分段線性配準結果
Figure10.
The result of piecewise linear registration
表1
分段線性配準前后圖譜與圖像的歸一化互信息值
Table1.
Normalized mutual information values of atlas and image before and after piecewise linear registration
對于本實驗所使用的四組數據，配準前的圖譜與圖像的歸一化互信息的均值為 1.060 7，而在配準后二者的歸一化互信息均值為 1.114 2，通過本章所提出的基于腦皮層標識點的分段仿射配準算法，實現了圖譜與圖像之間的初步配準，使二者之間的互信息增加了 5.04%，取得了較好的結果。對分段線性變換后的圖譜使用 B 樣條自由形變模型進行非剛性變換后的結果如圖 11 所示，使用更為精細的網格劃分方式進行對比，可見圖譜與患者圖像的內部紋理也更加趨于一致。
圖11
對分段線性配準所得圖像使用非剛性配準的結果
Figure11.
Deformable registration results for piecewise linear registration images
根據表 2 可知，分段線性配準后圖譜與圖像的歸一化互信息均值為 1.114 0，而在分段線性配準基礎上對二者繼續進行非剛性配準，所得的歸一化互信息值為 1.129 2，非剛性配準后二者的歸一化互信息值增長了約 1.36%。隨后，將圖譜中的配準結果應用于患者術前圖像以及術后圖像，即可實現如圖 12 所示的術前核團位置評估結果。分別計算術前核團位置評估結果以及專家在術前圖像序列中根據經驗值標記結果的重心，所得的物理坐標見表 3。
圖12
核團標記信息對比
Figure12.
Comparison of nuclear marker information
表3
算法評估結果與專家評估結果核團重心坐標
Table3.
Nucleus center of gravity coordinates of the algorithm evaluation results and expert evaluation results
表2
非剛性配準前后圖譜與圖像的歸一化互信息值
Table2.
Normalized mutual information values of atlas and image before and after deformable registration
經過計算，算法結果與專家根據經驗標記計算得出的重心誤差為（1.59 ± 1.063）mm。將術前圖像與術后圖像進行配準，并將配準所得的形變向量場應用至術前核團位置評估結果，可得到術后核團位置評估結果。隨后，對術后圖像中的電極進行手動分割，對術后核團與電極進行可視化，可得到如圖 13 所示的核團與電極三維可視化效果。
圖13
丘腦底核與電極術后三維可視化結果
Figure13.
Postoperative three-dimensional visualization results of nuclei and electrodes
3 總結
本文結合江蘇省人民醫院實際 DBS 影像數據特點，充分使用影像數據構建了 DBS 由術前到術后的效果評估管線，同時結合前人的研究成果，引入了人腦圖譜對核團位置進行估計，達到了較為標準的預測結果。整個配準流程都使用互信息作為相似性測度，適用于 MR 圖像不同的加權結果。在國內外現有的研究中，較多采用將患者手術前后的影像配準至圖譜中，以圖譜作為參考依據，會引入部分誤差，而本文為每例患者都構造了個性化的圖譜，得出的核團以及電極位置是基于患者本身坐標空間，相比于參考文獻[2]更多地利用了患者自身的信息，所得的結果對相應的患者更有代表性。在處理術后數據時，將大腦作為彈性模型，消除了手術前后大腦形變產生的誤差。
本文使用的數據在縱向的空間分辨率較低，因此也給結果帶來了一定的誤差。在后續研究中，需要考慮提升數據的質量。本文使用的互信息測度具有較好的魯棒性，但是有很大的運算量，由于本文所采用的數據為三維 MR 序列，因此該管線處理所需的運算時間較長，需要 3 個小時左右，如需將本文的方法用于術中實時評估，這是最需要解決的問題，也是后續研究的重要方向。
引言
腦深部電刺激（deep brain stimulation，DBS）手術，對于帕金森病、特發性震顫等疾病是一種有效的治療手段。在帕金森病治療中，DBS 主要通過在丘腦底核植入電極進行刺激，從而達到恢復患者運動功能的效果。患者在接受各項評估并達到手術要求后可以進行 DBS 手術治療。其第一步是進行手術計劃，患者將在局部麻醉條件下安裝手術頭架，并進行頭部磁共振掃描（magnetic resonance imaging，MRI）以及計算機斷層掃描（computed tomography，CT），通過影像融合確定立體定向頭架和電極植入點的定位。隨后是進行植入手術，根據手術計劃在患者頭皮標記點位置鉆圓孔，將探測電極按照計劃路徑植入腦內，術中通過刺激器測試電生理信號來進行微調。這一過程中，電極植入的精確度直接關系著手術的效果。最后是進行術后程控，在此過程中出現的失語、認知障礙以及疼痛等副作用也與電極植入精確度有關。因此，術后對于電極植入位置的準確評估是一個迫切需要解決的問題。
目前，臨床上一般使用高分辨 MRI 獲得患者術前的全腦信息以及核團位置。在手術后，由于電極為金屬材質，不能使用高分辨 MRI 對植入部位進行成像，因此一般使用 CT 或低分辨 MRI 對患者頭部進行掃描。在臨床中，通常采用 T2 加權方式對丘腦底核進行成像。由于丘腦底核體積較小，且與周圍組織沒有明顯的邊界，因此直接通過觀察判斷其位置會引入誤差。其次，術后的影像由于分辨率較低，只能在其中觀察到電極位置。同時，患者在接受術前術后掃描時無法保證體位的一致。這三方面的因素為準確判斷術后電極與核團的相對位置帶來了困難。在實際操作中，T2 加權成像方式需要的時間較長，往往只能得到核團區域的腦部圖像，這對后續判斷以及算法的設計提出了更高的要求。
DBS 術后效果評估的主要難點為患者丘腦底核的精確分割，人腦圖譜在核團分割過程中起到很大作用。Horn 等[1]將術前術后的影像數據同時向人腦圖譜配準，通過圖譜中的特定標記獲得丘腦底核的位置，并在此基礎上實現了 DBS 術后評估的工具包。張芷齊等[2]在標準腦圖譜空間下將雙側丘腦底核刺激治療患者的高分辨率術前圖像與低分辨率的術后圖像進行配準，并重建 DBS 電極的路徑，為醫生的術后程控提供輔助指導。Videen 等[3]根據術后 CT 圖像中的電極金屬偽影定位電極中心，并將手術前后圖像與 Mai 圖譜配準。上述方法在一定程度上可以解決丘腦底核分割以及電極定位問題，但是由于腦圖譜和患者腦部存在差異，因此核團分割結果的誤差難以預測。da Silva 等[4]將患者數據使用非剛性配準方法配準至標準圖譜，解決了患者與圖譜結構上的差異帶來的誤差，但是在配準過程中會導致電極的扭曲，會在結果中引入誤差。郝斌[5]對丘腦底核進行手動標記，實現了術后丘腦底核區域可視化，并從術后 MRI 數據中得到電極位置，分析電極植入位置與 DBS 手術療效之間的關系。本文使用非剛性配準方法，將腦圖譜轉化為與患者匹配度較高的個性化圖譜，減少了術后數據到圖譜的配準過程，從而減少在此過程中引入的誤差，同時根據術前術后的數據特點設計了數據處理流程，實現了一個端到端 DBS 術后效果評估管線，即一套將 DBS 手術圖像用于術后效果評估的標準流程。
本文旨在利用帕金森病 DBS 術前高分辨圖像中的核團信息及術后低分辨圖像中的電極信息，將二者進行配準，從而實現電極與核團空間位置的可視化，為術后程控提供可靠信息。同時，在研究中引入標準腦圖譜，將圖譜與術前圖像配準，解決術前影像中核團難以分割的問題。在實現術前影像中丘腦底核分割后，將術前與術后的圖像配準，從而得到電極與核團的相互關系。
1 算法設計
本文根據上述 DBS 手術影像數據的特點，在確保充分利用信息的情況下設計了如圖 1 所示的影像評估管線，首先將術前 T1 以及 T2 數據配準至同一空間，同時將術前 T1 數據以及配準后的 T2 數據標準化至 MNI 坐標空間，得到標準空間下的 T1 以及 T2 數據。隨后將標準的 ICBM152 圖譜與標準空間下的 T1 數據經過分段線性配準以及非剛性配準，得到患者的個性化圖譜。最后為了避免手術導致的大腦形變引入誤差，使用非剛性配準算法將患者術后 T1 數據與術前的個性化圖譜進行配準，根據圖譜中的標記確定術后核團的位置，并且將圖像中的電極提取出來，從而分析二者的相對位置關系。
圖1
影像評估管線流程圖
Figure1.
Flow chart of image evaluation pipeline
1.1 圖像預處理
首先，對患者的 MRI 圖像進行配準，將術前 T2 加權的 MRI 圖像以及術后 T1 加權的 MRI 圖像分別與術前 T1 加權的 MRI 圖像使用分割及配準工具包（Insight Segmentation and Registration Toolkit，ITK）[6]中的剛性配準算法進行配準，將手術前后的數據標準化到同一坐標空間中。在剛性配準過程中，由于組織的邊緣以及灰度在不同的加權方式下無法保持一致，因此采用互信息作為圖像的相似性度量。互信息是信息論中的一個基本概念，通常用于描述兩個系統間的統計相關性，可以用熵來表達[7]。
系統 A 的熵定義為
|
兩個系統的聯合熵為
|
兩個系統間的互信息為
|
當兩幅圖像空間位置完全一致時，可以用其中一幅圖像表達另一幅圖像的信息，可以使用聯合概率分布以及完全獨立的概率分布間的廣義距離估計互信息：
|
1.2 患者圖像與圖譜配準
在完成配準過程后，所有患者圖像都被標準化到同一坐標空間內。由于核團在患者 MR 影像中不明顯，需要通過將患者 T2 加權圖像與 MNI-ICBM152 圖譜配準實現對核團區域的識別。為了完成這一目標，首先需要將患者的腦圖像使用的坐標空間轉換成與圖譜中人腦方向一致的 MNI 坐標空間[8]。MNI 坐標空間以大腦前聯合（anterior commissure，AC）、后聯合（posterior commissure，PC）和腦中隔面為基準，與過 AC-PC 的連線的水平面（簡稱 AP 平面），以及過 AC 和 PC 且垂直于腦中隔和水平面的兩個平面（簡稱 VAC 平面及 VPC 平面）組成的平面將大腦分為 12 個長方體區域，如圖 2 所示。
圖2
MNI 坐標空間示意圖
Figure2.
MNI coordinate space diagram
完成坐標轉換需要識別腦中隔面，并將其旋轉至垂直位置，本文采用的是 Volkau 等[9]提出的使用 K-L 距離提取腦中隔面的方法。該方法是一種基于統計學的三維方法，魯棒性較好。K-L 距離衡量的是相同事件空間里的兩個概率分布的差異情況。設兩個隨機變量 p = {pi}和 q = {qi}，其中 pi 和 qi 表示第 i 種事件的發生概率。K-L 距離的定義為：
|
腦中隔是通過大腦半球間裂的平面，是含有腦脊液最多的矢狀面。如圖 3 所示，中隔面上腦脊液灰度值出現的概率最大，K-L 距離反映了兩個切片上灰度值出現概率的變化情況，由于中隔面的變化最大，因此與其他矢狀面的 K-L 距離最大。
圖3
距離中隔面兩側 2 cm、1 cm 的矢狀面與中隔面的比較圖
Figure3.
Comparison of the mid-sagittal plane and the two sides of the septum 2 cm, 1 cm from the plane
將中隔面校正至垂直后需要將 AP 平面校正至水平面，根據中隔面上的 AC、PC 點計算出需要圍繞旋轉的單位向量（nx，ny，nz）以及角度 θ，根據羅德里格斯旋轉公式，旋轉矩陣為：
|
由于人腦的個體差異較大，若直接與腦圖譜執行精確配準算法，容易導致算法在優化時陷入局部極值影響最終結果，因此需要采用分段仿射變換，將腦圖譜中的人腦大小變換成與患者較為接近的大小，為后續精確配準算法提供一個良好的初值。使用 Liu 等[10]提出的方法確定 AC、PC 點。根據 MNI 標準空間定義，如圖 4 所示，在 AP 平面上找出：過中隔面的額葉最前點 A；過中隔面的枕葉最后點 P；過 VPC 的左臨時葉最邊緣點 L；過 VPC 的右臨時葉最邊緣點 R。在 VAC 平面上找出過中隔面的腦組織最下點 I，在 VPC 上找出過中隔面的頂葉最上點 S。根據這些點計算出患者與圖譜根據 MNI 空間定義的 12 個分塊在每個維度上的線性比例系數[11]。根據所得系數對腦圖譜進行分段縮放，隨后計算縮放后圖譜與患者 MR 圖像中的 AC、PC 點位置的差值，對整個圖譜進行平移，使圖譜與患者前后聯合重合。
圖4
AP、VAC、VPC 上的標識點
Figure4.
The identification points on the AP, VAC and VPC
1.3 圖譜與患者圖像間的非剛性配準
在對圖譜完成上述變換后，模擬人腦與患者腦部的大小基本達到了一致，但還存在局部失配的現象，因此需要進行更為精細的非剛性配準。非剛性配準通過使標準圖譜形變擬合患者圖像，使差異較大的不同個體圖像間達到配準。本文所采用的是基于 B 樣條的自由形變模型（free form deform，FFD）的三維非剛性配準[12]。
基于 B 樣條函數的配準方法主要利用控制點和 B 樣條函數來描述非線性幾何變換。移動位于樣條函數上的控制點，使樣條函數彎曲，樣條函數之間的空間也會產生扭曲，達到自由形變的效果。在三維空間中基于 B 樣條函數的自由形變可以描述為三個一維三次 B 樣條函數的張量積：
|
其中 P 表示控制點網格，βi 表示 i 階 B 樣條函數，具體形式如下
|
由于圖譜與患者 MR 圖像間的配準屬于多模態圖像配準，因此同樣采用互信息作為相似性測度。
1.4 手術前后核團位置評估
標準的人腦圖譜根據以上變換轉化為與患者匹配度更高的個性化圖譜，ICBM-152 圖譜包含了如圖 5 所示的丘腦底核精確分割信息，因此將計算所得的形變網格應用至核團標識像中，便可獲得患者核團位置的評估結果。
圖5
丘腦底核在圖譜中的標記
Figure5.
Mask of the subthalamic nucleus in the atlas
術后患者影像雖然在圖像預處理時標準化到了術前患者圖像坐標空間中，但是由于電極植入過程需要打開患者顱腔，顱內壓力變化以及隨后的電極植入過程都會導致大腦形變[13]，因此需要對術前圖像向術后圖像進行非剛性配準。
2 實驗結果
2.1 實驗對象
受試者為 4 名帕金森患者，均已知曉并同意本研究使用其臨床數據。受試者在江蘇省人民醫院完成雙側丘腦底核 DBS 手術，植入電極為美國美敦力公司 3389 電極。手術前，使用德國西門子公司 3.0T MRI 系統，在快速自旋回波序列 T1 相重復時間（time of repetition，TR）1 600 ms、回波時間（time of echo，TE）2.48 ms、可視范圍（field of view，FOV）250 mm × 250 mm、像素 512 × 512 以及層厚 1.5 mm 條件下，以矢狀位無間隔掃描獲取全腦數據，隨后在快速自旋回波序列 T2 相 TR 5 000 ms、TE 84.00 ms、FOV 229 mm × 229 mm、像素 384 × 384、層厚 2 mm 條件下以橫斷位掃描獲取核團部位數據。手術后，使用美國通用電氣公司 1.5T MRI 系統，在快速自旋回波 T1 相 TR 4 000 ms、TE 105.62 ms、FOV 250 mm × 20 mm、像素 512 × 512、層厚 2 mm 條件下以橫斷位掃描獲得帶有電極的 MRI 圖像。
本文采用基于蒙特利爾神經科學研究所標準立體空間標準圖譜 MNI-ICBM152。該圖譜是將 152 名志愿者的高空間分辨率掃描結果，通過仿射變換與 MNI305 進行對應，再將 152 個大腦進行平均所得到的標準模板，已經被國際大腦映射聯盟（International Consortium for Brain Mapping，ICBM）作為標準模板。本文實驗算法設計及實驗過程在 Matlab 以及 Visual Studio 2010 上完成，并使用 ITK 中的部分算法。
2.2 原始數據
原始數據如圖 6 所示。術前數據為用于觀察大腦全局形態的矢狀位 T1 全腦數據和用于觀察核團的 T2 橫斷位數據（由于儀器條件所限，T2 加權方式掃描過程較長，因此只提供核團部位的數據）。為了防止金屬電極對患者造成傷害，術后數據只提供核團部位的橫斷位低分辨率數據。
圖6
原始臨床數據
Figure6.
The original clinical data
2.3 圖像預處理結果
剛性配準結果如圖 7 所示。可見算法執行后，術前術后核團部位與術前全腦數據完好地配準到了同一坐標空間中。
圖7
不同圖像與術前 T1 加權 MR 剛性配準結果
Figure7.
Rigid registration results of different images with preoperative T1-weighted MR
中隔面校正算法的執行結果如圖 8 所示，可見該算法較好地將中隔面校正到了與 x 軸垂直的狀態。
圖8
中隔面校正算法執行前（左側）后（右側）
Figure8.
Images before (left) and after (right) the middle septal correction
圖 9 為將圖像序列按照 AC-PC 點旋轉的前后對比效果，其中右側圖像中的患者大腦已校正至 MNI 標準空間。
圖9
AP 平面校正執行前（左側）后（右側）
Figure9.
Images before (left) and after (right) the AP plane correction
2.4 圖譜與患者圖像配準結果
由各組圖像的配準結果可以看出，ICBM152 圖譜的腦部較為狹長，將配準前后的圖譜與圖像進行融合，可得到如圖 10 所示的分塊融合結果。從中可以看到，經過配準的圖譜與臨床數據的腦組織在大小上更加一致。為了對分段仿射配準的效果進行量化，本文計算了配準前后圖譜與圖像之間的歸一化互信息值，其結果如表 1 所示。歸一化互信息值的區間為[1，2]，互信息值越大代表兩幅圖像間的配準效果越好。
圖10
分段線性配準結果
Figure10.
The result of piecewise linear registration
表1
分段線性配準前后圖譜與圖像的歸一化互信息值
Table1.
Normalized mutual information values of atlas and image before and after piecewise linear registration
對于本實驗所使用的四組數據，配準前的圖譜與圖像的歸一化互信息的均值為 1.060 7，而在配準后二者的歸一化互信息均值為 1.114 2，通過本章所提出的基于腦皮層標識點的分段仿射配準算法，實現了圖譜與圖像之間的初步配準，使二者之間的互信息增加了 5.04%，取得了較好的結果。對分段線性變換后的圖譜使用 B 樣條自由形變模型進行非剛性變換后的結果如圖 11 所示，使用更為精細的網格劃分方式進行對比，可見圖譜與患者圖像的內部紋理也更加趨于一致。
圖11
對分段線性配準所得圖像使用非剛性配準的結果
Figure11.
Deformable registration results for piecewise linear registration images
根據表 2 可知，分段線性配準后圖譜與圖像的歸一化互信息均值為 1.114 0，而在分段線性配準基礎上對二者繼續進行非剛性配準，所得的歸一化互信息值為 1.129 2，非剛性配準后二者的歸一化互信息值增長了約 1.36%。隨后，將圖譜中的配準結果應用于患者術前圖像以及術后圖像，即可實現如圖 12 所示的術前核團位置評估結果。分別計算術前核團位置評估結果以及專家在術前圖像序列中根據經驗值標記結果的重心，所得的物理坐標見表 3。
圖12
核團標記信息對比
Figure12.
Comparison of nuclear marker information
表3
算法評估結果與專家評估結果核團重心坐標
Table3.
Nucleus center of gravity coordinates of the algorithm evaluation results and expert evaluation results
表2
非剛性配準前后圖譜與圖像的歸一化互信息值
Table2.
Normalized mutual information values of atlas and image before and after deformable registration
經過計算，算法結果與專家根據經驗標記計算得出的重心誤差為（1.59 ± 1.063）mm。將術前圖像與術后圖像進行配準，并將配準所得的形變向量場應用至術前核團位置評估結果，可得到術后核團位置評估結果。隨后，對術后圖像中的電極進行手動分割，對術后核團與電極進行可視化，可得到如圖 13 所示的核團與電極三維可視化效果。
圖13
丘腦底核與電極術后三維可視化結果
Figure13.
Postoperative three-dimensional visualization results of nuclei and electrodes
3 總結
本文結合江蘇省人民醫院實際 DBS 影像數據特點，充分使用影像數據構建了 DBS 由術前到術后的效果評估管線，同時結合前人的研究成果，引入了人腦圖譜對核團位置進行估計，達到了較為標準的預測結果。整個配準流程都使用互信息作為相似性測度，適用于 MR 圖像不同的加權結果。在國內外現有的研究中，較多采用將患者手術前后的影像配準至圖譜中，以圖譜作為參考依據，會引入部分誤差，而本文為每例患者都構造了個性化的圖譜，得出的核團以及電極位置是基于患者本身坐標空間，相比于參考文獻[2]更多地利用了患者自身的信息，所得的結果對相應的患者更有代表性。在處理術后數據時，將大腦作為彈性模型，消除了手術前后大腦形變產生的誤差。
本文使用的數據在縱向的空間分辨率較低，因此也給結果帶來了一定的誤差。在后續研究中，需要考慮提升數據的質量。本文使用的互信息測度具有較好的魯棒性，但是有很大的運算量，由于本文所采用的數據為三維 MR 序列，因此該管線處理所需的運算時間較長，需要 3 個小時左右，如需將本文的方法用于術中實時評估，這是最需要解決的問題，也是后續研究的重要方向。
