張 斌，謝松云，李寧飛

（西北工業(yè)大學 電子信息學院，陜西 西安 710129）

EEG與MEG都能夠提供具有良好時間分辨率 （通常為毫秒級）的數(shù)據(jù)，但空間分辨率有限。與之相對應(yīng)，功能磁共振（function Magnetic Resonance Imaging， fMRI）提供了良好的空間分辨率卻相對較差的時間分辨率。研究人員們開始從將腦電與功能磁共振技術(shù)進行結(jié)合的角度尋求靈感和算法[1]，其目的就是結(jié)合兩者的優(yōu)勢，并最終得到更為合理精確的源定位結(jié)果。人腦建模是腦磁圖（Magnetoencephalogram，MEG）與腦電圖（Electroencephalogram，EEG）計算中正逆問題的先決條件。同時它也是基于EEG/fMRI的多模態(tài)醫(yī)學數(shù)據(jù)融合中約束方法首先要解決的問題。實際應(yīng)用中廣泛應(yīng)用的往往是球頭模型，通過將人腦部簡化為球使得正逆問題得到簡化。最為常用的為三層球頭模型，該模型包含頭皮、顱骨以及大腦三層，其它的類似球模型包括五層球模型以及重疊球模型等。然而，在使用偶極子方法進行求逆時，由于球模型與真實頭部的差異，在相關(guān)的數(shù)值運算以及激活源定位時會引入巨大的誤差[1]。這一缺陷也使得球模型的應(yīng)用受到了限制。為了提高新模型的準確程度，研究者開始考慮在構(gòu)建模型時結(jié)合人腦MRI數(shù)據(jù)。

隨后產(chǎn)生了幾種建立真實頭模型的數(shù)學方法，如基于體元的有限元方法（Finite Element Method，F(xiàn)EM）以及基于面元的邊界元方法（Boundary Element Method， BEM）[2]，這兩種方法都是通過使用三角元來構(gòu)建模型。為了保證模型的準確性，使用的三角元必須很小，因此需要大量的節(jié)點坐標，隨之而來的就是龐大的運算量，致使真實頭模型無法得到廣泛應(yīng)用[3]。本文提出了一種基于非均勻有理B樣條（Non-Uniform Rational B-Splines，NURBS）的構(gòu)建真實頭模型的新方法。NURBS曲面是一種在工業(yè)上早已得到廣泛應(yīng)用的方法，其它領(lǐng)域也有許多基于NURBS的建模方法研究。插值與擬合是最為基礎(chǔ)的兩種方法[4]。這里使用對MRI圖像進行插值的方法重建真實頭模型。該方法在保證結(jié)果良好的情況下建模速度更快且所使用的點更少。

1 總體流程設(shè)計

整個建模流程如圖1所示，首先對MRI數(shù)據(jù)進行圖像分割、邊緣提取等預(yù)處理，接著將圖像邊緣數(shù)據(jù)信息進行整合并規(guī)范化，最后對預(yù)處理得到的數(shù)據(jù)插值，反求出真實頭模型。

圖1 插值法構(gòu)建頭模型流程圖Fig.1 Flow chart of head modeling by interpolation

2 圖像預(yù)處理

2.1 圖像分割

圖像預(yù)處理部分的第一步是對MRI圖像進行分割[5]。使用的數(shù)據(jù)為Compressed NifTI（.nii.gz）格式，是磁共振常用的數(shù)據(jù)格式。掃描圖像為217×181×181矢狀圖，即從如圖2所示的矢狀面角度看，每層圖像大小為217×181，共掃描了181層。

由于實際圖像不可能在特定結(jié)構(gòu)區(qū)域有穩(wěn)定的灰度值，因此需要應(yīng)用一定的形態(tài)學方法對圖像加以處理，去除可能的噪點，得到邊界后才能進行下一步的初始化。

圖2 圖像分割Fig.2 Image segmentation

2.2 邊緣提取

該模塊的第二部是邊界提取，這里使用變形方法完成頭部MRI圖像邊緣提取。該方法首先對原始MRI圖像進行閾值處理并二值化如圖3（a）所示。二值化后求得該圖像的重心和半徑，利用其重心和半徑構(gòu)造一個位于頭內(nèi)部的初始球面，再通過不斷變形迭代更新曲面，最終得到所要求的邊緣結(jié)果。接下來再進行形態(tài)學處理，使用圖像閉合方法，即對圖像進行先膨脹后腐蝕，再進行填孔運算，最終得到圖3（b）所示的大腦邊界圖像

圖3 邊界提取Fig.3 Boundary extractioin

3 基于NURBS的曲面重構(gòu)

3.1 數(shù)據(jù)初始化

經(jīng)過圖像的預(yù)處理模塊，我們已經(jīng)得到了MRI圖像的頭部邊界，要以這些邊界點作為已知數(shù)據(jù)點進行曲面插值。要進行曲面插值，首先要求這些數(shù)據(jù)點必須是有序的，所以在對MRI圖像進行分割得到邊緣點后，必須按照一定規(guī)則，通過邊界跟蹤的方法按順序排列這些邊緣點，從而為插值做準備。此外，由于邊緣提取方法所限，不可避免的會出現(xiàn)一些孤立的噪聲點，可能導(dǎo)致跟蹤方法失敗，因此在邊緣跟蹤開始前，需要先對提取的邊緣數(shù)據(jù)去噪，即將孤立點以及與邊緣相連的一些異常點刪除。

邊界跟蹤即從圖像的一個邊緣點出發(fā)，按照一定的規(guī)則搜索下一邊緣點直到滿足特定終止條件，由此得到目標邊界。具體步驟為：首先，確定起始搜索點。本文取圖像左上角第一個有值點作為起始搜索點S，將其位置存入邊界點序列并作為第一個邊界點。將所有點搜索狀態(tài)置0，存儲在一個與圖像等大小的零矩陣中，且該點不是已搜索得到的邊界點。

然后，通過對提取的腦邊緣進行八鄰域搜索，將各掃描層大腦邊界點的坐標數(shù)據(jù)按順序以矩陣形式存儲在一個元胞數(shù)組內(nèi)。由于總的數(shù)據(jù)點數(shù)目依然很大，如果直接將這些點用于構(gòu)建模型仍然會給計算機很大的運算負載，與此同時，生成的模型也缺乏光順性。因此需要在進行曲面重構(gòu)前對數(shù)據(jù)進行預(yù)處理。

由于MRI數(shù)據(jù)中各層腦切面大小不同，因此所存儲的邊界數(shù)據(jù)點也不相同，首先必須將數(shù)據(jù)矩陣規(guī)范為相同大小，并在這一過程中降低矩陣中元素的數(shù)目。本文所用MRI數(shù)據(jù)中共有146層矢狀面包含大腦數(shù)據(jù)，故邊界矩陣中包含146行數(shù)據(jù)點。各行對應(yīng)數(shù)據(jù)由數(shù)十列到數(shù)百列不等。如果進行曲面重構(gòu)時需要一個大小為m×n的矩陣，那么對于數(shù)據(jù)長度超過n列的行來說就可以進行采樣，而對于少于n列數(shù)據(jù)點的則按順序插入一些點。通過這個方法就可以生成一個m×n的邊界點矩陣，表示為 Qij（i=0，1，…，m；j=0，1，…，n）。 由于最終建模時所需要的n小于多數(shù)切面對應(yīng)的邊界數(shù)據(jù)點，因此在對數(shù)據(jù)進行預(yù)處理后的總數(shù)據(jù)量有大幅減少。這里n的選取可以根據(jù)各層所需要的重構(gòu)曲線質(zhì)量自由選擇。

3.2 曲面重構(gòu)

接下來的問題就是如何根據(jù)邊界數(shù)據(jù)使用NURBS重構(gòu)大腦曲面。由于數(shù)據(jù)點已經(jīng)被規(guī)范化為一個矩陣，因此這個問題可以通過一般的曲面插值得到解決。

一個（p，q）次的B樣條曲面表示形式如下[6]：

其中 Pi，j為控制點或定義點，Ni，q（u）與 Nj，q（v）為 p 與 q 次下的B樣條基函數(shù)，u和v為節(jié)點，對應(yīng)節(jié)點向量分別為：

曲面插值問題可以表示為

其中

這里使用雙三次樣條曲面插值，即p=3，q=3。插值過程包含以下幾步：

1）數(shù)據(jù)點參數(shù)化，采用積累弦長參數(shù)化方法

首先令d為總弦長

那么有u0=0 un=1

該方法可用于大多數(shù)情況，且對應(yīng)產(chǎn)生的插值曲線結(jié)果光順性較好，符合通常的要求，因此是應(yīng)用最為廣泛的參數(shù)化方法。它在數(shù)據(jù)點分布不均勻或控制多邊形弦長不等長的情況下，能夠避免如上均勻參數(shù)化方法的問題，給出符合數(shù)據(jù)點或弦長分布的參數(shù)化結(jié)果，所以它在某種程度上也被粗略的認為是弧長參數(shù)化。這里需要強調(diào)的是，曲線插值最后生成的結(jié)果光順性好壞并不僅僅取決于所選的參數(shù)化方法，同時還由所選的插值方法所決定。此外，這里采用弦長參數(shù)化方法也并不等同于說最終的插值曲線參數(shù)是積累弦長。

2）計算節(jié)點向量U

節(jié)點通?？梢赃M行等距離分割

然而通常并不推薦使用這種方法，這是因為當結(jié)合積累弦長參數(shù)化方法使用時，有可能會產(chǎn)生一個奇異方程組。因此使用如下的平均方法

通過使用這一方法可以使得節(jié)點的分布反映出參數(shù)化后點uk的分布。

3）利用公式（4）對U向使用NURBS曲線插值，求得一組“中間”控制點。

4）根據(jù)公式（5）對V向使用NURBS曲線插值，求得控制點 Pi，j。

5）使用 U，V 以及 Pi，j計算反求 NURBS 曲面。

4 重構(gòu)結(jié)果

整個建模過程在MATLAB環(huán)境下完成。最終重構(gòu)真實頭模型中的大腦部分如圖4所示。

圖4 大腦重構(gòu)結(jié)果Fig.4 Reconstructed result

計算機仿真結(jié)果如圖4所示，最終的生成結(jié)果良好。除此之外，在數(shù)據(jù)預(yù)處理時，數(shù)據(jù)點的個數(shù)由超過70 000點降至7 000余點。而如果使用三角元構(gòu)建模型，數(shù)字將大幅提升。上述建模過程在一臺處理器為Intel Core（TM）2 Duo Processor，內(nèi)存為2G的普通臺式機上僅耗時十秒左右。

5 結(jié) 論

該建模系統(tǒng)使用模塊化思想結(jié)合NURBS新方法構(gòu)建真實腦模型，具有高速、少數(shù)據(jù)建模的特點。對于該方法有效性的更進一步評估可以通過在腦電的正逆問題計算中進行量化分析。此外，未來開展的研究中可以使用更為復(fù)雜的NURBS方法。NURBS作為真實頭模型構(gòu)建中一種很有潛力的方法，是值得投入更多精力的。

[1]Lalancette M，Quraan M，Cheyne D.Evaluation of multiplesphere head models for MEG source localization[J].Physics in Medicine and Biology，2011，56（17）:5622-5636.

[2]Akahn Z，Acar C E，Gencer N G.Development of realistic head models for electromagnetic source imaging of the human brain engineering[J].Proceedings of the 23rd Annual International Conference of the IEEE Medicine and Biology Society，2001，1:899-902.

[3]He J J，Shen H，Hu D W.A survey of the EEG/fMRI fusion analysis:head models，methods and applications[J].Computer Engineering&Science，2007，29（12）:74-81.

[4]Piegl L，Tiller W.The NURBS Book[M].New York:Springer，1997.

[5]Smith S M.Fast robust automated brain extraction[J].Human Brain Mapping，2002，17（3）:143-155.

[6]Piegl L，Tiller W.Reducing control points in surface interpolation[J].IEEE Computer Graphics and Applications，2000，20（5）：70-74.

[7]李焱，時芝勇，海曉濤.基于改進蟻群算法的配電網(wǎng)重構(gòu)[J].陜西電力，2010（9）：22-25.LI Yan,SHI Zhi-yong,HAI Xiao-tao.Distribution Network Reconfiguration Based on Improved Ant Colony Algorithm[J].Shaanxi Electric Power，2010（9）：22-25.