張 斌，謝松云，李寧飛

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

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

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

1 總體流程設(shè)計(jì)

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

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

2 圖像預(yù)處理

2.1 圖像分割

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

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

圖2 圖像分割Fig.2 Image segmentation

2.2 邊緣提取

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

圖3 邊界提取Fig.3 Boundary extractioin

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

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

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

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

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

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

3.2 曲面重構(gòu)

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

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

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

曲面插值問題可以表示為

其中

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

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

首先令d為總弦長(zhǎng)

那么有u0=0 un=1

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

2）計(jì)算節(jié)點(diǎn)向量U

節(jié)點(diǎn)通?？梢赃M(jìn)行等距離分割

然而通常并不推薦使用這種方法，這是因?yàn)楫?dāng)結(jié)合積累弦長(zhǎng)參數(shù)化方法使用時(shí)，有可能會(huì)產(chǎn)生一個(gè)奇異方程組。因此使用如下的平均方法

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

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

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

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

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

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

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

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

5 結(jié) 論

該建模系統(tǒng)使用模塊化思想結(jié)合NURBS新方法構(gòu)建真實(shí)腦模型，具有高速、少數(shù)據(jù)建模的特點(diǎn)。對(duì)于該方法有效性的更進(jìn)一步評(píng)估可以通過在腦電的正逆問題計(jì)算中進(jìn)行量化分析。此外，未來開展的研究中可以使用更為復(fù)雜的NURBS方法。NURBS作為真實(shí)頭模型構(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]李焱，時(shí)芝勇，海曉濤.基于改進(jìn)蟻群算法的配電網(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.