基于超像素的多模態(tài)MRI腦膠質(zhì)瘤分割

蘇坡, 楊建華, 薛忠

(1.西北工業(yè)大學 自動化學院, 陜西 西安 710129; 2.Houston Methodist研究所, 美國 休斯頓 77030)

腦膠質(zhì)瘤(glioblastoma multiforme, GBM)是腦腫瘤中最常見的、死亡率最高的一種惡性腫瘤。統(tǒng)計表明，40%的腦腫瘤患者為腦膠質(zhì)瘤。腦膠質(zhì)瘤患者術(shù)后存活時間的中位數(shù)僅為8個月，而5年以上的存活率幾乎為零[1]。GBM在多模態(tài)MRI圖像中呈現(xiàn)出一片異質(zhì)的腫瘤區(qū)域。這片區(qū)域通常包括3個部分：壞死腫瘤(necrosis)、活動腫瘤(enhanced tumor)以及腫瘤擠壓周圍正常腦組織所形成的水腫(edema)。由于GBM腫瘤在組織形態(tài)上的復雜性與特殊性，單模態(tài)MRI無法清晰地反映GBM的不同組織結(jié)構(gòu)。相比之下，多模態(tài)MRI圖像中含有豐富的組織結(jié)構(gòu)信息，被廣泛用于GBM的診斷和治療。本文的多模態(tài)MRI主要包括T2(T2加權(quán)成像)、T1PRE(T1加權(quán)成像)、T1POST(T1增強成像)和FLAIR (液體衰減反轉(zhuǎn)成像)。在不同的模態(tài)下，GBM組織圖像呈現(xiàn)出不同的特征：活動腫瘤在T1POST呈現(xiàn)高信號，壞死部分在T1POST呈現(xiàn)低信號，而水腫部分在T2和FLAIR呈現(xiàn)高信號。GBM的分割是指根據(jù)這些特征，把GBM組織從正常的腦組織中標記或分割出來。

文獻中，基于像素或者體素的分割算法廣泛應(yīng)用于GBM的分割。這類算法的基本思想是根據(jù)每個像素在多模態(tài)圖像上亮度信息、紋理信息等把該像素點分類到相應(yīng)的類別中。分類的算法包括無監(jiān)督的聚類和有監(jiān)督的學習。例如，Clark等[2]提出了一種基于FCM(fuzzy C-means algorithm)的模糊均值聚類的算法。該算法以多模態(tài)圖像的灰度作為特征向量，首先利用FCM對所有體素點進行聚類得到初始的分類，然后根據(jù)對稱性，灰度分布等先驗知識對初始分類進行優(yōu)化，得到最終的分割結(jié)果。由于FCM聚類時，沒有考慮空間鄰域信息，并且GBM組織的灰度分布會產(chǎn)生重疊，因此容易產(chǎn)生誤分割。

基于圖形分割的算法[3-4]現(xiàn)在也非常流行。這類算法用圖的頂點來描述圖像的像素，用圖的邊描述2個像素的相似性，由此形成一個網(wǎng)絡(luò)圖，然后通過解決能量最小化問題把圖分割成子網(wǎng)絡(luò)圖，使不同子網(wǎng)絡(luò)圖之間的差異和同一子網(wǎng)絡(luò)圖內(nèi)部的相似性達到最大。這類算法通常需要解決一個求解廣義特征向量問題，當圖像比較大時，這類算法會遭遇計算復雜度大的問題。除了以上2類算法，基于Level set[5]的分割算法也廣泛應(yīng)用于GBM分割，但是由于GBM組織灰度不均勻，并且GBM組織之間經(jīng)常沒有明顯的邊界，采用這類算法容易出現(xiàn)邊緣泄露的問題。

最近，超像素或超體素(supervoxel)[6-7]引起了很大的關(guān)注，基于超像素和超體素的分割算法也成為研究熱點。相比于傳統(tǒng)的像素和體素，超像素和超體素具有計算效率高，符合人類視覺感知，能高效表示圖像信息等優(yōu)點。

為了充分利用超像素的優(yōu)點，提高GBM分割的精確性和魯棒性，本文提出了基于超像素的多模態(tài)MRI GBM分割算法。首先通過帶加權(quán)距離的局部k-均值聚類算法，把多模態(tài)MRI圖像過分割成一系列均勻、緊湊并能很好地貼合圖像邊緣的超像素。然后，利用動態(tài)區(qū)域合并算法對超像素進行逐步合并。最后通過后處理得到GBM的最終分割結(jié)果。實驗表明本文算法能取得較好的分割效果。

1 研究方法

1.1 算法框架

圖1給出本文的算法框架。算法包含4步:①圖像預處理，包括多模態(tài)圖像的共同配準和去腦殼；②利用本文提出的帶加權(quán)距離的局部k-均值聚類算法把多模態(tài)圖像分割為一系列亮度均勻并能很好地貼合圖像邊緣的超像素；③利用動態(tài)區(qū)域合并算法對②產(chǎn)生的超像素進行合并。通過區(qū)域合并生成若干塊均勻的有意義的區(qū)域(合并后的不同區(qū)域用不同的顏色表示，見圖1)；④根據(jù)亮度分布等對合并后的區(qū)域進行后處理，最終完成對GBM的分割。

圖1 算法框架

1.2 預處理

預處理包括多模態(tài)圖像的共同配準(co-registration)和去腦殼。圖像的共同配準保證多模態(tài)圖像同一位置的像素對應(yīng)腦部同一組織。本文利用FSL[8]的FLIRT工具完成對多模態(tài)MRI圖像的共同配準。選取T2作為參考圖像，把其他模態(tài)的圖像配準到T2上。配準采用剛性變換并以互信息量作為圖像相似性測度。

去腦殼是腦圖像處理的一個常見步驟。它一方面可以減少后續(xù)處理的計算量，另一方面也會減少非腦實質(zhì)組織對后續(xù)處理的影響。本文利用一個模板圖像和它去腦殼后的圖像，通過配準的方法實現(xiàn)去腦殼操作。首先把該模板腦圖像通過仿射變換配準到T2上。接著用第一步配準產(chǎn)生的仿射變換矩陣對去腦殼后的模板圖像進行變換就實現(xiàn)了對T2的去腦殼操作。

1.3 基于局部k-均值聚類的超像素生成

在文獻[6]中，Achanta等提出了一種基于局部k-均值聚類算法，把彩色圖像過分割為一系列均勻、緊湊并能很好貼合圖像邊緣的超像素。該算法的基本思想是根據(jù)像素點的圖像特征和空間位置特征定義一個新的距離函數(shù)，然后根據(jù)新距離函數(shù)對像素點進行局部的k-均值聚類從而產(chǎn)生超像素。該算法主要輸入?yún)?shù)為K，即待生成的超像素的個數(shù)。首先，在待分割的彩色圖像均勻分布K個初始聚類中心點。根據(jù)K和整個圖像像素點總數(shù)目N可以計算出初始聚類中心點之間的間距S：

(1)

然后，選取像素在CIELAB顏色空間的顏色值[lab]T和自身空間位置[xy]T組成5維特征向量[labxy]T。與傳統(tǒng)的k-均值聚類采用歐式距離作為相似性測度不同，Achanta定義了一個新的相似性測度，該相似性測度考慮了所產(chǎn)生超像素的大小。新的相似性測度定義像素點i到第k個聚類中心的距離為

(2)

式中

(3)

(4)

(liaibixiyi)T表示第i個像素的特征向量，i∈[1,N]。(lkakbkxkyk)T表示第k個聚類中心的特征向量，k∈[1,K]。df為像素點i到第k個聚類中心的圖像特征距離，測量2個像素亮度特征的相似性。dxy為像素i到第k個聚類中心的空間距離，測量2個像素空間位置的相似性。S為公式(1)中初始聚類中心的之間的空間間隔，它的引入相當于引入了超像素大小的信息。參數(shù)m控制所產(chǎn)生超像素的緊湊程度。m越大，空間距離dxy占主導，所產(chǎn)生的超像素就越緊湊；反之，m越小，亮度特征距離df占據(jù)主導，產(chǎn)生的超像素能更好地貼合圖像邊緣，但其大小和形狀不是很緊湊。傳統(tǒng)的k-均值聚類搜索范圍為整個圖像空間，為了減少k-均值聚類的計算復雜度，Achanta限定搜索范圍為聚類中心周圍2S×2S區(qū)域。這就是該算法被稱為局部k-均值聚類的原因。

本文對Achanta的算法進行了擴展，實現(xiàn)從多模態(tài)MRI產(chǎn)生超像素。選取像素在T2、T1PRE、FLAIR和T1POST的亮度值[f1f2f3f4]T和空間位置[xy]T組成6維特征向量[f1f2f3f4xy]T。為了產(chǎn)生更好的超像素，我們對Achanta算法進行改進，提出了一種帶加權(quán)距離的局部k-均值聚類。我們的依據(jù)是不同模態(tài)的圖像對GBM分割的貢獻不同，如T2和FLAIR含有比較多的水腫信息，T1POST含有較多的GBM壞死部分和活動部分信息。對T1POST模態(tài)進行強調(diào)，定義新的圖像特征距離df：

(5)

參數(shù)ω(0<ω<1)為T1POST模態(tài)的權(quán)重，ω=0.25時各個模態(tài)圖像的權(quán)重相等(標準局部k-均值聚類)，ω>0.25時T1POST模態(tài)占圖像特征距離的比重大于其他模態(tài)。ω的引入可以實現(xiàn)對超像素進行細調(diào)，使超像素在保持緊湊性的同時還能更好地貼合GBM活動腫瘤和壞死部分的邊緣?？臻g距離dxy的選取跟Achanta的算法一致。

圖2 采用標準局部k-均值聚類和帶加權(quán)距離的局部k-均值聚類產(chǎn)生超像素的對比圖

圖2給出采用標準局部k-均值聚類(等權(quán)重距離)和加權(quán)距離k-均值聚類產(chǎn)生超像素的例子。其中設(shè)置K=300,m=70。從圖2的箭頭所指的超像素可以看出采用帶加權(quán)距離k-均值聚類算法產(chǎn)生的超像素能更好地貼合GBM壞死部分和活動腫瘤部分的邊緣。

1.4 動態(tài)區(qū)域合并算法

上一節(jié)中，多模態(tài)圖像被過分割成一系列均勻、緊湊并能很好貼合圖像邊緣的超像素。如圖3所示，多模態(tài)MRI被過分割成一系列超像素，參數(shù)設(shè)置K=600,m=70,ω=0.4。圖3第一排圖像依次為T2、T1PRE、FLAIR和T1POST模態(tài)圖像；第二排和第三排顯示用本文算法產(chǎn)生的超像素。

圖3 多模態(tài)MRI圖像被過分割成超像素

我們采用區(qū)域合并的方法對超像素進行處理進而實現(xiàn)對GBM的分割。為了方便起見，用區(qū)域鄰接圖(region adjacency graph, RAG)模型描述所產(chǎn)生的超像素。令G(V,E)表示一個無向圖，vi∈V(i=1,2,…,M)為圖的一個頂點，代表第i個超像素。E?V×V表示相鄰的2個頂點的邊，邊的權(quán)重w(vi,vj)表示相鄰的2個超像素i和j之間的相似度。

區(qū)域合并的過程就是根據(jù)區(qū)域鄰接圖頂點的相似性對頂點進行合并同時對相應(yīng)的數(shù)據(jù)進行更新的過程。區(qū)域合并算法成功的關(guān)鍵在于能否解決2個問題：①設(shè)計區(qū)域合并的規(guī)則，即2個區(qū)域滿足什么樣的條件就對其合并；②設(shè)計合并停止條件，即區(qū)域合并到什么程度就停止區(qū)域合并。針對這2個問題，研究者提出了不同的區(qū)域合并算法。Nock等[10]采用相似性統(tǒng)計檢驗方法實現(xiàn)區(qū)域合并：若2個相鄰的區(qū)域具有統(tǒng)計意義上的相似性，則對這2個區(qū)域進行合并，反之不進行合并。然而現(xiàn)有的大部分區(qū)域合并算法均存在容易產(chǎn)生過合并，欠合并或者兩者混合的情況。針對現(xiàn)有區(qū)域合并算法的不足，Peng等[10]提出了一種動態(tài)區(qū)域合并算法并證明了該算法產(chǎn)生的區(qū)域合并結(jié)果能保持一定的全局屬性。該算法把區(qū)域合并看做是一個由相似性度量和一致性度量2個部分組成的推斷問題。相似性度量解決尋找候選合并區(qū)域的問題，一致性度量則決定是否對這些候選區(qū)域進行合并。在一致性度量部分，Peng引入序貫概率比檢驗(sequential probability ratio test, SPRT)方法，通過SPRT來確定候選區(qū)域是否具備一致性。實驗結(jié)果表明Peng的算法優(yōu)于其他區(qū)域合并算法。本文采用Peng的動態(tài)區(qū)域合并算法來實現(xiàn)區(qū)域合并。定義區(qū)域合并謂詞P為:

P(v1,v2)=

v1和v2表示RAG中2個相鄰的超像素v1和v2。w(vi,vj)為RAG 2個頂點vi和vj的相似度。Ω1和Ω2分別表示與v1和v2相鄰的超像素的集合。謂詞P(v1,v2)表示當且僅當同時滿足條件(a)：v1和v2互為最相似的鄰居；條件(b)：v1和v2具有一致性時對v1和v2進行合并，否則不進行合并。條件(b)暗含停止區(qū)域合并的條件。若沒有條件(b)，所有的超像素最終將合并成為一個區(qū)域。定義2個超像素vi和vj相似度w(vi,vj)為：

(6)

N(vi)和N(vj)分別表示超像素vi和vj中所含像素的個數(shù)。vi和vj均為向量，分別表示超像素vi和vj在4種模態(tài)圖像的亮度平均值向量。

對于條件(b)，Peng采用序貫概率比檢驗(SPRT)的方法來進行判斷。定義零假設(shè)H0∶θ=θ0，即2個相鄰的超像素vi和vj具有一致性；備擇假設(shè)H1∶θ=θ1，即vi和vj不具有一致性。根據(jù)假設(shè)檢驗第一類錯誤的概率α和第二類錯誤的概率β確定2個閾值A(chǔ)和B(B<0

(7)

推斷問題定義的斷言是隱狀態(tài)，不能直接觀測，但它在統(tǒng)計上與可以直接觀測到的信息x相關(guān),這里x為多維隨機變量代表超像素中的像素在各個模態(tài)圖像的亮度。定義序貫似然比和序貫似然比累加和分別為：

(8)

式中:P0(xi|θ0)和P1(xi|θ1)為第i次抽樣觀測到x的條件概率，這里用高斯分布模型來近似x的分布，這2個條件概率可以表示成：

(9)

Ia和Ib均為向量，分別表示對待合并的相鄰超像素vi和vj進行隨機抽樣(抽樣個數(shù)為各自所含像素數(shù)目的一半)后，各自樣本4種模態(tài)亮度的均值；Ia+b為這2個超像素所組成的大區(qū)域中所有像素4種模態(tài)亮度的均值；CI表示這2個超像素所組成的大區(qū)域所有像素4種模態(tài)亮度的協(xié)方差矩陣；λ1和λ2為2個標量參數(shù)。

為了解決SPRT檢驗中抽樣檢驗次數(shù)的不確定性問題，采用截斷型SPRT檢驗算法，即預設(shè)一個最壞情況抽樣次數(shù)的上限N。然后根據(jù)下面的規(guī)則對假設(shè)的2種模式進行判斷：

其中n表示第n次抽樣。本文所有實驗SPRT參數(shù)均設(shè)為α=0.05,β=0.05,N0=10,λ1=10,λ2=1。

根據(jù)區(qū)域合并的條件(a)和條件(b)反復地對相鄰的超像素進行合并和更新直至所有兩兩相鄰的超像素不能合并為止，就完成了整個區(qū)域合并過程。

1.5 后處理

通過動態(tài)區(qū)域合并，超像素最終被合并為大大小小的幾十個區(qū)域。由于區(qū)域合并是對相鄰的區(qū)域進行合并，而實際情況是GBM同一組織空間上可能是不相鄰的，因此需要對合并結(jié)果進行后處理進而得到最終的GBM分割。

首先根據(jù)亮度和空間位置剔除掉與背景相鄰的大腦外圍的脊髓液區(qū)域和去腦殼沒有去除干凈的非腦實質(zhì)組織區(qū)域。脊髓液在FLAIR模態(tài)為低亮度，而未被去除干凈的非腦實質(zhì)組織的區(qū)域所含的像素數(shù)目比較少。因此可以事先設(shè)定FLAIR亮度閾值和像素數(shù)目閾值剔除以上組織。接著用區(qū)域中所有像素在4種模態(tài)的平均亮度作為區(qū)域的4種模態(tài)亮度，見圖4a)。我們對區(qū)域的4種模態(tài)亮度進行Kernel PCA[13]降維，把四維降至二維。圖4b)為對區(qū)域4種模態(tài)亮度進行Kernel PCA降維后的結(jié)果，其中實心圓點代表水腫區(qū)域，實心方塊代表活動腫瘤部分，實心三角代表壞死部分。然后對區(qū)域亮度Kernel PCA降維后的結(jié)果進行k-均值聚類，再根據(jù)不同組織在不同模態(tài)圖像亮度的分布規(guī)律就可以得到最終的分割結(jié)果見圖5。k-均值聚類時，我們可以從Kernel PCA的結(jié)果觀察得到k-均值聚類的類別數(shù)k，例如在根據(jù)圖4b)設(shè)定聚類數(shù)目k=7、8。注意在聚類的過程中，同時加入如下約束：若k-均值聚類時2個相鄰的區(qū)域被分為一類，則不把這2個區(qū)域劃為一類，因為動態(tài)區(qū)域合并的結(jié)果表明它們不具有一致性。如圖4b)所示，距離2個實心圓點最近的那個空心點，按照標準k-均值聚類，這個空心點極有可能跟這2個實心圓點歸為一類，但是由于它和其中的一個實心圓點在實際空間位置上相鄰的，根據(jù)附加約束這個點不能和實心圓點歸為一類。

圖4 區(qū)域合并后的區(qū)域亮度進行Kernel PCA結(jié)果

圖5 最終分割結(jié)果(灰度值從高到低依次為水腫、壞死和活動腫瘤)

2 實驗結(jié)果與分析

為了評價新算法的性能，本文對15個GBM病人的多模態(tài)MRI腦圖像數(shù)據(jù)進行定性和定量實驗。每個病人多模態(tài)MRI圖像包括T2，T1PRE，F(xiàn)LAIR和T1POST，每種模態(tài)的圖像大小均為256×256×21, 圖像的分辨率為0.78 mm×0.78 mm×6.5 mm。本文算法的性能依賴于超像素的質(zhì)量。超像素主要由K，m和ω控制。通過反復實驗，我們設(shè)定K=600,m=70,ω=0.4。

為了定量評估新算法的性能，我們把本文算法和其他2種經(jīng)典分割算法進行比較，這2種經(jīng)典算法是：基于FCM算法[2]和Ncut[3]。

我們選用Dice相似性系數(shù)評價分割算法性能。Dice相似性系數(shù)的定義為：

(10)

A表示醫(yī)生手動分割的GBM組織像素的集合，B表示分割算法分割出來的GBM組織像素的集合。Dice相似性系數(shù)越大，表明分割的結(jié)果越好。圖6給出本文算法和基于FCM算法，Ncut算法的性能比較圖。從圖6可以看出和其他2種經(jīng)典算法相比，本文的算法更加精確。

圖6 算法性能比較圖

3 結(jié) 論

本文提出了一種基于超像素的多模態(tài)MRI腦膠質(zhì)瘤分割算法：把超像素和動態(tài)區(qū)域合并這2種算法結(jié)合到一起，充分利用2種算法各自優(yōu)點，提高GBM分割的魯棒性和準確性。對15個GBM病人數(shù)據(jù)的對比實驗證明了本文算法的先進性。

參考文獻：

[1] Patel M R, Tse V. Diagnosis and Staging of Brain Tumors[J]. Seminars in Roentgenology, 2004, 39(3): 347-360

[2] Clark M C, Hall L O, Goldgof D B, et al. Automatic Tumor Segmentation Using Knowledge-Based Techniques[J]. IEEE Trans on Medical Imaging, 1998, 17(2): 187-201

[3] Shi J B, Malik J. Normalized Cuts and Image Segmentation[J]. IEEE Trans on Pattern Analysis and Machine Intelligence, 2000, 22(8): 888-905

[4] Corso J J, Sharon E, Dube S, et al. Efficient Multilevel Brain Tumor Segmentation with Integrated Bayesian Model Classification[J]. IEEE Trans on Medical Imaging, 2008, 27(5): 629-640

[5] Ho S, Bullitt E, Guido G. Level-Set Evolution with Region Competition: Automatic 3-D Segmentation of Brain Tumors[C]∥IEEE International Conference on Pattern Recognition, 2002: 532-535

[6] Achanta R, Shaji A, Smith K, et al. SLIC Superpixel Compared to State-of-the-Art Superpixel Methods[J]. IEEE Trans on Pattern Analysis and Machine Intelligence, 2012, 34(11): 2274-2282

[7] Lucchi A, Smith K, Achanta R, et al. Supervoxel-Based Segmentation of Mitochondria in EM Image Stacks with Learned Shape Features[J]. IEEE Trans on Medical Imaging, 2012, 31(2), 474-486

[8] Smith S M, Jenkinson M, Woolrich M W, et al. Advances in Functional and Structural MR Image Analysis and Implementation as FSL[J]. Neuroimage, 2004, 23: 208-219

[9] Nock R, Nielsen F. Statistical Region Merging[J]. IEEE Trans on Pattern Analysis and Machine Intelligence, 2004, 26(11): 1452-1458

[10] Peng B, Zhang L, Zhang D. Automatic Image Segmentation by Dynamic Region Merging[J]. IEEE Trans on Medical Imaging, 2011, 20(12): 3592-3605

[11] Sch?lkopf B, Smola A, Müller K R. Nonlinear Component Analysis as a Kernel Eigenvalue Problem[J]. Neural Computation, 1998,10(5):1299-1319