基于Laplace方程的腦皮層厚度測量方法

邰瀅瀅，梁永霞

(遼寧大學 信息學院，遼寧 沈陽 110036)

0 引言

大腦是人體最重要、最復雜的器官，它控制著人體各組織器官功能的正常運作.大腦皮層是大腦最外層的薄皮狀結構，是人腦最高級的中樞，因此很多疾病的產生也會通過大腦皮層腦部形態學特征變化體現出來.通過對大腦皮層形態特征的變化規律進行研究，可以高效地診斷和治療疾病，甚至進一步預防疾病的發生，這對于現代醫學發展有著十分重要的意義.人眼視力的好壞也是與大腦皮層結構有著密不可分的關系的.視覺剝奪會導致大腦皮層形態結構的變化，研究發現這些變化既包括局部腦結構的改變，如皮層厚度、灰質體積等，也包括腦結構網絡的改變.針對這種情況，如果能夠發現近視者早期特征和腦部形態特征變化規律，不僅對視力下降者早期治療非常重要，也對預防近視起著決定性的作用.先天盲人和正常人在腦皮層形態學上都存在著某些差別.通過研究這兩類人群的腦部差異，能夠有助于盲人的康復治療，甚至能夠使盲人通過其他的直覺補償，來得到更好的生活體驗.因此，對大腦皮層進行形態學分析有助于發現疾病病理特征和發病規律.大腦皮層形態學指標有很多種，其中腦皮層厚度就是形態學研究中的一個重要指標.通過估計整個大腦皮層的厚度信息，可以發現視力剝奪引起的大腦皮層的形態學變化.

1 相關知識介紹

腦皮層厚度是指灰質皮層外表面到白質皮層外表面之間的距離，皮層厚度的測量方法主要分為基于網格表面測量和基于體數據測量兩大類[1].在基于體數據的測量方法中有由Rosenfeld和Pfaltzc[2]提出的基于距離變換的方法.這種方法在進行距離變換時是在二維圖像切片上實現的，所以在某些程度上會使得測量結果產生偏差.Lohmann等[3]考慮了二維圖像距離變換時會丟失部分三維信息，導致大腦組織結構的特征變化，其算法在三維結構上進行，保證了拓撲不變性，使結果更合理.Jones等[4]提出了一種新的皮層厚度定義方法，對腦皮層構建勢能場模型，并且用拉普拉斯方程來描述這個勢能場，通過求解這個偏微分方程，獲取各個等勢面，再根據等勢面計算出場線.場線位于皮層內表面和外表面之間，與所有經過的等勢面都垂直，這樣我們就可以把皮層內過任一點的場線的長度定義為該點的厚度.盡管這種求解二階偏微分方程的方法計算量很大，但這種方法是目前較為科學的一種測量方法.

基于網格表面的測量方法中有一種比較常見的方法由Fischl等[5]在2000年提出.該方法具有良好的魯棒性.其原理是基于大腦皮層曲面的網格頂點，計算了一次外表面(灰質)到內表面(白質)頂點的最短距離，一次內表面(白質)到外表面(灰質)網格頂點的最短距離，然后求兩次最短距離的平均值作為該點的皮層厚度.這種測量方法的基礎在于從分割的MRI圖像重建出三維網格皮層表面.因此這類方法測得的結果與網格重建和圖像分割等過程緊密相關.例如求圖1中A點處的皮層厚度，首先通過三維重建得到網格化的腦皮層表面，再基于網格頂點計算每一點的皮層厚度.具體步驟是先計算A到B的最短歐氏距離，再計算B到C的最短歐氏距離，取AB+BC的平均值就得到A點的皮層厚度值.

圖1 腦皮層厚度算法示意圖

2 基于網格的腦皮層厚度測量

基于網格的腦皮層厚度測量的基礎是從MRI二維圖像上重建出網格狀的三維皮層表面，再由Free Surfer直接計算出網格上每一點的皮層厚度.這其中需要經過一系列的預處理步驟，如圖2所示.

圖2 圖像預處理過程圖

2.1 圖像預處理

圖像處理包括預處理和迭代計算兩部分.在進行分析之前要經過一系列的預處理，圖像的預處理包含很多過程，對于不同的后處理其對應的預處理過程也不相同.總的來說，圖像預處理大致包括以下幾步：去噪、配準、不均勻場校正、組織分割、內外皮層曲面重建等.這些預處理過程都是在醫學圖像處理中非常熱門的研究方向.在本文中主要是大腦皮層的形態學指標的計算和統計，因此，預處理的質量決定著形態學指標測量的準確度.圖像的預處理步驟如下：

2.1.1 去噪

圖像去噪是指減少數字圖像中噪聲的過程.現實中的數字圖像在數字化和傳輸過程中常受到成像設備與外部環境噪聲干擾等影響，稱為含噪圖像或噪聲圖像.圖像中噪聲的來源有許多種，這些噪聲來源于圖像采集、傳輸、壓縮等各個方面.在這里的去噪主要是去除核磁共振成像過程中產生的噪聲信號和人腦其他組織信號，這些信號對后期進行大腦皮層的處理沒有作用，甚至會造成誤差影響，因此需要去除掉這些噪聲信號.

2.1.2 配準

配準是一個類似于標準化的處理，將個體的核磁共振圖像通過平移、旋轉、縮放使其與標準模板對齊，并最大化與模板相似的一個過程.通過配準，可以得到一個與模板圖像一一對應的個體圖像，通過模板又可以得到不同個體圖像間的對應關系，從而使得不同空間內的圖像之間相互比較成為可能.

2.1.3 不均勻場校正

在核磁共振成像的過程中，磁場強度不能完全保證均勻一致，造成了圖像灰度值的不均勻性.同一腦組織(灰質、白質、腦脊液)的灰度值和對比度在不同區域相差很大.這對于主要依靠灰度值或對比度來進行分割算法的正確率有很大的影響，進而直接影響到后續的處理和分析.因此，必須對不均勻場進行校正.

2.1.4 腦組織分割

腦組織的分割是預處理中必要的步驟，對于后續計算和分析有重要作用.腦組織分割一般首先把大腦分割成指定的幾類，如灰質、白質和腦脊液三類.首先根據圖像的灰度值進行簡單的預劃分，將體素按照灰度值進行分類.選取出那些最可能屬于某一類的體素.對于無法從灰度上區分的剩余體素，再利用大腦幾何信息來進行精細分類.

2.1.5 曲面重建

在完成腦組織分割之后，我們就可以根據分割結果對大腦灰質和白質外邊界進行三維曲面重建.曲面重建就是要恢復出隱藏在三維體數據中的腦皮層結構的幾何信息.因為從一組二維斷層圖像上理解病人身體內部器官的三維結構是非常困難的，即使是經驗豐富的醫生也很難準確地描述器官的立體結構和確定病變體的方位.另外，三維曲面重建技術的出現也給醫學教學領域帶來了很大的幫助，學生可以通過三維重建獲得直觀立體的人體器官和組織圖像，形象地理解和掌握各個器官和組織的構成與功能以及常見病變的發生情況，提高學習的質量和效率.

3 基于體數據的腦皮層厚度測量

基于體特征的腦皮層厚度估計方法是在三維立方體體素網格上進行運算的，并通過求解拉普拉斯方程來得到腦皮層厚度[6].這種方法在計算前也是需要進行預處理的.因為大腦是三維結構的，直接對二維圖像進行分析會丟失部分三維信息，所以要預處理轉化為三維體素結構來進一步測量.圖3為基于體數據的皮層厚度測量流程圖.在基于體數據的測量中，本文的核心思想是在盡可能準確分割的基礎上，借助二階偏微分方程測量大腦皮層厚度.

圖3 數據處理流程圖

3.1 圖像分割

由于大腦形態結構極其復雜，僅使用一種分割方法來分割圖像并不能得到較好的分割結果，因此本文作者將多種分割方法結合起來，發揮各自的優點，中和單一分割方法帶來的不足.首先根據圖像的灰度值進行簡單的預劃分[9]，將體素按照灰度值進行分類.選取出那些最可能屬于某一類的體素.對于無法從灰度上區分的剩余體素，再利用大腦幾何信息來進行精細分類.本文測量方法的流程圖如圖3所示.

3.2 腦皮層厚度測量

通過一系列預處理提取出腦皮層區域，并定義出皮層的內外邊界，根據數學原理用拉普拉斯方程法測量厚度，主要包括以下幾步：在腦皮層邊界建立拉普拉斯方程得到勢能場；計算單位向量場；計算皮層厚度.傳統體數據厚度的測量采用拉普拉斯方程求解方法，在整個三維腦皮層空間內進行迭代求解，迭代公式為

(1)

ψi+1(x，y，z)是第i次迭代時坐標(x，y，z)處像素的勢能值.當求解得到ψ之后，可以由式(2)求解得到單位向量場線T，其中?為梯度算子.

(2)

然后，就可以用偏微分方程計算皮層厚度[11].假設皮層區域內每個點的單位向量場線都是由內表面S指向外表面S′，則-T表示沿同一條場線從外表面S′指向內表面S.這里定義兩個長度函數L0和L1，L0表示從內表面到某點的場線長度，L1表示從外表面到某點的場線長度.

根據微分幾何，L0和L1都滿足下面的一階線性偏微分方程

(3)

那么任一點x處的腦皮層厚度值即為穿過該點的場線長度L0和L1的和：

L(x)=L0(x)+L1(x)

(4)

其中L(x)就是點x處的皮層厚度值.

4 實驗與分析

4.1 實驗環境

本實驗運行的計算機配置為2.4 GHz Intel Core i3和4 GB RAM.通過上面的一系列預處理步驟，得到一個基于三角網格的大腦皮層表面，實驗用FreeSurfer軟件來進行MRI二維斷層圖像的預處理.FreeSurfer軟件是用于MRI圖像處理和分析的一款開源軟件，可以進行多種模態數據的預處理、皮層重建等(http：//surfer.nmr.mgh.harvard.edu/)，該軟件只能在Linux系統下運行，所以選擇在Windows系統下安裝虛擬機，在虛擬機中安裝Ubuntu，然后將下載好的FreeSurfer安裝包解壓并安裝.然后計算出網格上的腦皮層頂點的厚度.本文中用于實驗的研究對象有先天性盲人組20例，年齡范圍為32～45歲，正常組20例，年齡范圍為35～49歲，均通過核磁共振掃描得到.選取的兩組對象的年齡、性別、智力狀態都相差不大，也沒有其他的疾病.部分相關信息如表1所示.對這40張二維圖像進行預處理得到三維數據.本文在分割步驟中使用SPM5軟件進行兩次粗分割和細分割，實現對所有實驗樣本的腦組織進行分割，然后在提取出的大腦皮層上建立拉普拉斯方程，依據3.2節的理論公式用迭代法求解方程，得到腦皮層厚度值.

表1 早期盲人樣本資料

(5)

4.2 實驗結果分析

經過對兩組人群腦皮層厚度的統計發現，早期盲人的視覺區皮層厚度總體上要比正常人視覺區的厚，如圖4所示，基于Laplace方程得到的皮層厚度統計，橫坐標是大腦皮層不同部位的學術名稱，縱坐標是皮層厚度值.從圖中可以看出，視力剝奪對于直接處理視覺信息的枕葉腦區影響尤為顯著，而且研究發現整個腦區在形態學結構上都與正常組有所不同，這也說明盲人在視覺區受到視力剝奪，為了彌補視力的不足就會在其他功能區域上得到強化[13].

圖4 大腦視覺區腦皮層統計圖

同時，本文分別使用兩種不同的計算方法得到正常組和早期盲人組的皮層厚度信息，表2是兩種方法在腦皮層不同部位測得厚度的均值，用同樣的檢驗方法進行了比較，得到對應的差異值q，結果顯示，在正常組和盲人組之間的差異中，基于網格的方法對應的q值為0.026 5，基于體數據的方法對應的q值為0.031 8.可以發現在基于體數據的方法中盲人和正常人之間的差異更大，這也說明本文準確分割的體數據方法有更高的測量精度.

表2 兩種方法在不同部位的皮層厚度均值

5 結束語

本文比較了基于網格和基于體數據的兩種計算腦皮層厚度的方法，并對早期盲人和正常人的皮層厚度進行了統計，結果表明基于體數據的測量方法具有更高的精度，但本文實驗所采集的被試樣本類別較少，后續可以依據不同的階段來進行分組，如先天盲、后天盲和正常人，可以對兩種方法進行多個對照組的比較，探究出盲人在先天和后天上腦皮層形態的變化.此外，關于腦部結構度量有很多形態學指標，后續我們可以計算更多數據指標，如皮層復雜度、灰質體積、灰質密度、白質體積等，通過分析這些指標可以更好地研究盲人的大腦結構，這對于醫學研究將具有十分重要的意義.