曹 敏，侯金秀，張?jiān)路?，鄧紅霞，李海芳

(太原理工大學(xué) 信息與計(jì)算機(jī)學(xué)院，山西 晉中 030600)

相關(guān)研究表明獼猴和人類的大腦結(jié)構(gòu)有著一定的相似性，獼猴(具有較高智力的非人靈長類動物)代表了最佳的用于驗(yàn)證假設(shè)的可侵入性模型[1]，相對于人腦，其在研究腦結(jié)構(gòu)和認(rèn)知方面有著無可比擬的優(yōu)勢。由于上述原因，近幾年關(guān)于獼猴的研究正在逐漸增多[2-3]。

對核磁圖像進(jìn)行分析是研究獼猴大腦的重要手段，作為核磁圖像分割領(lǐng)域的重要研究方向，核團(tuán)分割的準(zhǔn)確性及魯棒性在促進(jìn)腦研究領(lǐng)域的發(fā)展方面具有重要的研究價(jià)值和意義。然而對獼猴大腦的分割結(jié)果目前大多由研究人員借助人腦分割工具(如FSL)中提供的方法獲得，然后再進(jìn)行手動修正[4]。2017年，BALBASTRE[5]在高斯混合模型(gaussian mixture model，GMM)方法基礎(chǔ)上進(jìn)行改進(jìn)，應(yīng)用于獼猴的T2圖像進(jìn)行分割，但其能同時(shí)分割的腦組織有限。多圖譜分割方法與該方法相比可以引進(jìn)更多的先驗(yàn)信息，提高分割精度，所以很多研究人員將研究視線投向了多圖譜方法。一些研究人員發(fā)布獼猴的各種MRI模板、圖譜以及數(shù)據(jù)集[6-8]。MILHAM et al[9]對現(xiàn)有的非人靈長類公開數(shù)據(jù)集進(jìn)行了比較。陳偉導(dǎo)等[10]用多圖譜分割方法對獼猴腦組織進(jìn)行了分割。最近，隨著深度學(xué)習(xí)方法的迅速發(fā)展，有研究人員將深度學(xué)習(xí)技術(shù)應(yīng)用于核磁圖像分割。HUO et al[11]采用3D-UNET網(wǎng)絡(luò)對人腦進(jìn)行三維分割，ZHAO et al[12]利用深度學(xué)習(xí)對非人靈長類動物的大腦進(jìn)行提取。但此方法所需的訓(xùn)練數(shù)據(jù)集較大，還有較大的改進(jìn)空間。在多圖譜方法人腦分割領(lǐng)域，ZHANG et al[13]將稀疏表示的算法引進(jìn)到多圖譜標(biāo)簽融合的領(lǐng)域后，取得了較好的效果。嚴(yán)盟等[14-15]提出的方法在稀疏表示標(biāo)簽融合方面做出了一些努力。夏瑞等[16]提出了一種基于重采樣的改進(jìn)的多圖譜分割算法。LIU et al[17]提出了一種基于圖集配準(zhǔn)和線性化的核稀疏代表分類器的有效分割方法，但算法運(yùn)行時(shí)間較長。

本文選用稀疏表示多圖譜方法對獼猴腦組織的MR圖像進(jìn)行分割。為了提高大腦皮下核團(tuán)的分割準(zhǔn)確度，針對稀疏塊包含的信息較少和稀疏表示方法會丟失信息的問題，本文提出了一種改進(jìn)的稀疏表示多圖譜分割標(biāo)簽融合算法(multi-atlas segmentation label fusion algorithm based on improved sparse representation，ISRLF).該方法在稀疏表示的過程中引入圖譜標(biāo)簽，圖像中的標(biāo)簽信息用來增加圖像塊的先驗(yàn)信息，使計(jì)算的稀疏系數(shù)更加合理；然后通過改變信息熵的計(jì)算方式改進(jìn)互信息，將其用來衡量線性配準(zhǔn)后的圖譜與目標(biāo)圖像之間的相似度，將相似度的值與稀疏系數(shù)相結(jié)合計(jì)算出最終分配給圖譜的權(quán)重，然后進(jìn)行標(biāo)簽融合；隨后通過設(shè)計(jì)的一個(gè)相似性度量指標(biāo)將非局部塊加權(quán)融合的結(jié)果和稀疏表示融合的結(jié)果進(jìn)行綜合，該指標(biāo)包括Dice系數(shù)和余弦距離兩部分，最終達(dá)到了較好的分割效果。

1 背景描述

稀疏表示最初應(yīng)用于信號處理領(lǐng)域，其思想是將原信號用一些基本信號的組合代替，這樣可以盡可能地減少處理的數(shù)據(jù)量，如公式(1)所示，T為被表示信號，D為字典矩陣，α為稀疏系數(shù)。稀疏表示在圖像處理領(lǐng)域已被應(yīng)用于圖像的壓縮，超分辨率重建，特征提取表示等方面。稀疏系數(shù)α可利用公式(2)求得。式中x為待表示信號。

T=Dα.

(1)

(2)

2 ISRLF方法

2.1 方法框架

本文將稀疏表示多圖譜分割方法進(jìn)行了改進(jìn)，提出了一種將稀疏表示融合方法與非局部塊加權(quán)方法綜合在一起的新標(biāo)簽融合算法，對獼猴大腦皮下核團(tuán)進(jìn)行分割。圖1為本方法的框架圖。

算法步驟如下：

1) 對圖譜灰度圖像進(jìn)行線性配準(zhǔn)，得到線性配準(zhǔn)后的圖像Li.

2) 對線性配準(zhǔn)后的圖譜灰度圖像進(jìn)行非線性配準(zhǔn)，得到配準(zhǔn)后的結(jié)果Fi.將圖譜標(biāo)簽圖像非線性配準(zhǔn)到目標(biāo)圖像中，得到配準(zhǔn)后的結(jié)果Ii.

3) 對目標(biāo)圖像采用局部加權(quán)融合的方法進(jìn)行預(yù)分割，得到預(yù)分割結(jié)果。

4) 構(gòu)建目標(biāo)圖像標(biāo)簽塊Tpx和圖譜圖像標(biāo)簽塊Tpi，所有圖譜圖像標(biāo)簽塊cdi組成過完備字典D，即D=[Tp1,Tpi,Tpn]，進(jìn)行稀疏表示。

5) 引入互信息并對信息熵的計(jì)算方法進(jìn)行改進(jìn)，然后對目標(biāo)體素和線性配準(zhǔn)后的圖譜相對應(yīng)的平面進(jìn)行全局的相似性衡量，求平均值，改進(jìn)稀疏表示計(jì)算出的權(quán)重。

6) 分別用稀疏表示的標(biāo)簽融合方法和非局部塊加權(quán)的標(biāo)簽融合方法對目標(biāo)圖像進(jìn)行分割，得到各自的分割結(jié)果。

7) 采用Dice系數(shù)和余弦距離相結(jié)合的相似性度量標(biāo)準(zhǔn)對分割結(jié)果進(jìn)行度量，將兩種融合方法的分割結(jié)果再次進(jìn)行融合，最終得到分割結(jié)果。

圖1 方法框架圖Fig.1 Method framework

2.2 目標(biāo)圖像標(biāo)簽塊和圖譜圖像標(biāo)簽塊的構(gòu)建

一般的稀疏表示方法其圖像塊所包含的信息只有目標(biāo)圖像對應(yīng)位置的灰度值，沒有充分利用更多的先驗(yàn)信息。本方法將灰度圖像中的灰度信息和標(biāo)簽圖像中的標(biāo)簽信息結(jié)合，分別引入到目標(biāo)圖像和圖譜灰度圖像中構(gòu)成目標(biāo)圖像標(biāo)簽塊和圖譜圖像標(biāo)簽塊，使圖像塊攜帶的信息更加豐富，從而優(yōu)化標(biāo)簽融合時(shí)各圖譜的權(quán)重。

本文采用D99-SL、NMT和inia19三個(gè)圖譜進(jìn)行分割。為了將圖譜中的標(biāo)簽信息引入，需要對目標(biāo)圖像進(jìn)行預(yù)分割，本方法采用局部加權(quán)法對目標(biāo)圖像進(jìn)行預(yù)分割。構(gòu)建初始的僅含有灰度信息的圖像塊后，將對應(yīng)的標(biāo)簽值塊變成二值化的形式，再與初始圖像塊結(jié)合構(gòu)成新的圖像標(biāo)簽塊。二值化的具體方式為以中心體素的標(biāo)簽值為基準(zhǔn)，與之標(biāo)簽值一致的標(biāo)簽記為1，否則記為0.

圖譜圖像標(biāo)簽塊構(gòu)建的具體過程如圖2所示，以非線性配準(zhǔn)后的圖譜灰度圖像Fi中位于x處的體素作為中心，以r×r×r為大小構(gòu)建初始圖像塊adi，再從非線性配準(zhǔn)后的圖譜標(biāo)簽圖像It中，以對應(yīng)位置的標(biāo)簽為中心，以r×r×r為大小構(gòu)建標(biāo)簽塊bdi并二值化。將adi與bdi分別變成列向量后進(jìn)行拼接，組成圖譜圖像標(biāo)簽塊Tpi.

圖2 圖譜圖像標(biāo)簽塊的構(gòu)建Fig.2 Construction of atlas image label block

目標(biāo)圖像標(biāo)簽塊構(gòu)建的過程與圖譜圖像標(biāo)簽塊的構(gòu)建過程相似。首先，以目標(biāo)體素x為中心，r×r×r為大小構(gòu)建目標(biāo)圖像塊adx，然后以預(yù)分割結(jié)果的對應(yīng)位置x為中心，r×r×r為大小構(gòu)建標(biāo)簽塊bdx，分別變成列向量，最后將adx與bdx合并，組成目標(biāo)圖像標(biāo)簽塊Tpx.

對于目標(biāo)圖像標(biāo)簽塊Tpx，選取各圖譜圖像標(biāo)簽塊Tpi組成過完備字典D，即D=[Tp1,Tpi,Tpn].按照公式(2)進(jìn)行稀疏矩陣的計(jì)算，式中的x為目標(biāo)圖像標(biāo)簽塊，計(jì)算出的稀疏系數(shù)作為圖譜的初始權(quán)重矩陣W1.

2.3 基于互信息的圖像相似度計(jì)算

稀疏表示的標(biāo)簽融合方法只考慮了以目標(biāo)體素為中心的圖像塊中的體素，并沒有從全局的角度考慮圖譜圖像與目標(biāo)圖像的相似度，在標(biāo)簽融合時(shí)對整體圖像相似性的衡量有助于提高分割效果?；バ畔⑹切畔⒄撝械囊粋€(gè)概念，在圖像處理領(lǐng)域常通過衡量兩幅圖像的信息熵比較兩幅圖像的相似程度，所以本方法引入互信息度量參與融合的圖譜灰度圖像的相似度。它的值越大代表兩張圖片的相似性越高。公式(3)為圖像的信息熵計(jì)算公式，式中的概率pM由公式(4)求得，公式(4)中通過灰度直方圖進(jìn)行信息熵的計(jì)算，hi表示灰度值為i的體素?cái)?shù)量，分母表示該圖像的體素個(gè)數(shù)。

(3)

(4)

通?；叶葓D像之間的互信息是通過灰度直方圖進(jìn)行計(jì)算，對于核磁圖像來說有不足之處。例如，核磁圖像中灰度值為零的點(diǎn)表示的是圖像的背景，一般的互信息會將其納入相似度計(jì)算之中，而核磁圖像中只關(guān)注大腦區(qū)域，背景基本上可以忽略不計(jì)。假如核磁圖像中灰度值為零的點(diǎn)過多，而不為零的點(diǎn)過少，就可能會出現(xiàn)計(jì)算出的互信息值大，但真正重要的腦部結(jié)構(gòu)可能相似度并不高的情況。此外，只對灰度圖像進(jìn)行相似度計(jì)算引入的先驗(yàn)信息有限，不能很好地衡量兩幅圖像的相似性。

為了對上述兩種情況進(jìn)行改進(jìn)，本方法對信息熵的公式進(jìn)行了修改。首先將灰度值為零的體素剔除，不納入概率計(jì)算中，按照公式(5)計(jì)算灰度值的概率pM′.

(5)

其次由于大腦結(jié)構(gòu)中灰白質(zhì)占有著重要的位置，本方法關(guān)注的是皮下核團(tuán)，所以本方法通過配準(zhǔn)后的組織概率模板，按照公式(6)求得整體圖像中體素為灰質(zhì)的概率與體素為灰白質(zhì)概率的比值，達(dá)到衡量腦結(jié)構(gòu)體素相似度的目的，該計(jì)算結(jié)果作為參數(shù)β，引入到灰度圖像的信息熵計(jì)算之中。

(6)

為了增加先驗(yàn)信息，可以將標(biāo)簽圖像引入，將標(biāo)簽圖像引進(jìn)的方法是對標(biāo)簽圖像按照核團(tuán)的標(biāo)簽值進(jìn)行信息熵計(jì)算。標(biāo)簽圖像信息熵的計(jì)算方式如公式(7).Im表示標(biāo)簽圖像，m表示的是圖像中核團(tuán)的數(shù)量。由于Dice系數(shù)可以表征相同核團(tuán)的相似度，所以將其作為參數(shù)引入到計(jì)算過程中以增加先驗(yàn)信息。如公式(8)，求得目標(biāo)圖像分割結(jié)果IA與IB關(guān)于標(biāo)簽j的Dcie系數(shù)，將公式(8)計(jì)算出的值作為參數(shù)γ引入到改進(jìn)的信息熵的計(jì)算中。

(7)

γ=Dice(IAj,IBj) .

(8)

由此，經(jīng)過改進(jìn)的信息熵表達(dá)形式為公式(9).對于給出的兩幅圖像A和B，H1(A)和H1(B)分別為圖像的信息熵，計(jì)算方式為公式(10)、(11)，H1(A，B)為圖像A和B的聯(lián)合信息熵，計(jì)算公式為(12).I1(A,B)為A與B的互信息，計(jì)算公式為(13).

H1(M)=βH(M)+γH(IM) .

(9)

(10)

(11)

(12)

I1(A,B)=H1(A)+H1(B)-H1(A,B) .

(13)

將計(jì)算出的矢狀面、冠狀面、橫斷面三個(gè)不同面的互信息作為S1,S2,S3矩陣，提取矩陣中與目標(biāo)體素(a,b,c)對應(yīng)的相似度值，然后按照公式(14)進(jìn)行計(jì)算，得到體素x的相似度值，再按照公式(15)、(16)，得到最終的圖譜融合權(quán)重。

(14)

W2=W1×S(x) .

(15)

(16)

2.4 標(biāo)簽融合

基于圖像塊的融合方法有兩種，一種是通過一個(gè)相似性度量函數(shù)度量所有參與標(biāo)簽融合的圖像塊與目標(biāo)圖像塊的相似性的非局部塊加權(quán)融合方法，一種是稀疏表示方法，將目標(biāo)圖像塊由參與融合的圖譜圖像塊稀疏表示。

非局部塊加權(quán)的標(biāo)簽融合方法有著參與融合的圖像塊沒有被舍棄、冗余信息多的優(yōu)點(diǎn)，稀疏表示的融合方法與非局部塊加權(quán)的融合方法相比，該方法會舍棄一些信息。但該融合方法也有著非局部塊加權(quán)的融合方法難以比擬的優(yōu)勢，其在邊界點(diǎn)的分割效果上比非局部塊加權(quán)的融合方法有著更好的表現(xiàn)。

為了結(jié)合這兩種方法各自的優(yōu)點(diǎn)，提高對大腦皮下核團(tuán)的分割準(zhǔn)確度，本方法將非局部塊加權(quán)方法和稀疏表示的標(biāo)簽融合方法進(jìn)行了結(jié)合，從而提出了一種新的標(biāo)簽融合方法ISRLF.

基于圖像塊加權(quán)融合方法的過程如下：

1) 將各圖譜灰度圖像和標(biāo)簽圖像配準(zhǔn)到目標(biāo)圖像T中，得到配準(zhǔn)后的灰度圖像Fi和標(biāo)簽圖像Ii.

2) 提取T中以位置x為中心大小為r×r×r的目標(biāo)圖像塊Tpx.(本實(shí)驗(yàn)r值為3)然后進(jìn)行圖像塊的提取。提取過程為各圖譜以當(dāng)前要融合的體素x為中心，r×r×r作為搜索鄰域，提取所有的參與標(biāo)簽融合的圖像塊Tpi.

3) 采用某種方法計(jì)算出各圖像塊的權(quán)重，然后按照公式(17)進(jìn)行標(biāo)簽值概率值的計(jì)算。

4) 按照公式(18)得到標(biāo)簽的結(jié)果。

在本方法中提取的圖像塊為圖像標(biāo)簽塊，普通的非局部塊加權(quán)方法采用常用的歸一化相關(guān)系數(shù)函數(shù)(normalized correlation coefficient，NCC)衡量目標(biāo)圖像塊與提取出的圖譜圖像塊之間的相似度。稀疏表示的標(biāo)簽融合方法是構(gòu)建目標(biāo)圖像標(biāo)簽塊和圖譜圖像標(biāo)簽塊，將圖譜圖像標(biāo)簽塊作為過完備字典D=[Tp1,Tpi,Tpn]，通過公式(2)進(jìn)行稀疏系數(shù)的計(jì)算，然后通過公式(14)、(15)、(16)得到最終的各圖譜權(quán)重。然后按照公式(17)、(18)得到分割結(jié)果。

(17)

(18)

為了將兩種分割方法進(jìn)行融合，本方法設(shè)計(jì)了一個(gè)相似性度量函數(shù)。該函數(shù)包含兩個(gè)部分，第一部分為Dice系數(shù)，第二部分為余弦距離。公式中相似性度量結(jié)果與Dice系數(shù)成反比，與余弦距離成正比。記目標(biāo)體素x處預(yù)分割結(jié)果為R，非局部塊加權(quán)方法的融合結(jié)果為R1，稀疏表示融合方法的融合結(jié)果為R2.度量規(guī)則表示如下：

(19)

式中：Dice(R，Ri)為非局部塊加權(quán)(或者稀疏表示)融合方法的分割結(jié)果和預(yù)分割結(jié)果的Dice系數(shù)，d為常數(shù)；添加常數(shù)d是為了防止分母為零的情況出現(xiàn)；pi(x)表示非局部塊加權(quán)的結(jié)果或者稀疏表示融合方法的結(jié)果以體素x為中心的一個(gè)圖像塊轉(zhuǎn)化成的列向量，cos(p(x),pi(x))為兩個(gè)列向量的余弦距離；μ和1-μ代表兩者的權(quán)重；DRi值越大，說明相似度越低。

公式(20)為兩種方法的融合公式。當(dāng)非局部塊加權(quán)方法的分割結(jié)果和稀疏表示的分割結(jié)果的標(biāo)簽值相等時(shí)，就將該標(biāo)簽值賦予體素x，當(dāng)兩種方法的分割結(jié)果不相同時(shí)，通過公式(19)，計(jì)算相似度結(jié)果。若DR1>DR2，說明R1與預(yù)分割結(jié)果R之間的差距比R2與R之間的差距大，將R2的分割結(jié)果賦予體素x，若DR1

(20)

3 實(shí)驗(yàn)與結(jié)果分析

3.1 數(shù)據(jù)介紹

實(shí)驗(yàn)選用的測試數(shù)據(jù)集為牛津大學(xué)發(fā)布的獼猴數(shù)據(jù)集。該數(shù)據(jù)集用3T的掃描儀采集了20只雄性獼猴的數(shù)據(jù)。數(shù)據(jù)集包括T1，靜息態(tài)的fMRI。本實(shí)驗(yàn)使用T1數(shù)據(jù)。獼猴年齡分布為2.41歲～6.72歲(平均年齡為4.01歲)；重量分布為：4.35 kg～11.7 kg(平均重量6.57 kg)；體素分辨率為0.5 mm×0.5 mm×0.5 mm，TE為4.01 ms，TR為2 500 ms，TI為1 100 ms，翻轉(zhuǎn)角度為8°.

3.2 評價(jià)指標(biāo)

實(shí)驗(yàn)選用Dice系數(shù)評價(jià)分割質(zhì)量，其計(jì)算方法如公式(21)所示。公式中的A表示目標(biāo)圖像的金標(biāo)準(zhǔn)(即專家分割結(jié)果)，B表示用分割方法得出的分割結(jié)果。‖表示各集合的體素個(gè)數(shù)，|A∩B|表示集合A與集合B對應(yīng)體素標(biāo)簽相同的體素個(gè)數(shù)，|A|與|B|分別表示集合A與集合B中的體素個(gè)數(shù)。Dice系數(shù)取值范圍為0到1，越接近1，說明準(zhǔn)確率越高。

(21)

3.3 結(jié)果分析

本實(shí)驗(yàn)選取三個(gè)方法分別對海馬體、紋狀體、帶狀體這三個(gè)核團(tuán)的分割結(jié)果與金標(biāo)準(zhǔn)進(jìn)行比較。參與對比的方法有多數(shù)投票法(majority voting，MV)，聯(lián)合標(biāo)簽融合法(LW joint)以及本文的方法(ISRLF).

圖3(a)-(d)分別為目標(biāo)圖像和三個(gè)不同方法的分割結(jié)果。其中，3(a)為金標(biāo)準(zhǔn)分割結(jié)果，3(b)為MV方法的分割結(jié)果，3(c)為LW joint方法的分割結(jié)果，3(d)為ISRLF方法的分割結(jié)果。從圖中可以看出，LW joint方法的分割結(jié)果最差，ISRLF方法的效果最好。

圖3 三種不同方法的分割結(jié)果Fig.3 Segmentation results of three different methods

表1顯示了在牛津大學(xué)獼猴數(shù)據(jù)集上采用上述三種方法對選定核團(tuán)的分割準(zhǔn)確率。圖4為3個(gè)方法分割結(jié)果Dice系數(shù)的盒圖。從表1與圖4可以看出，MV方法雖然對于海馬體的分割準(zhǔn)確率比LW joint的準(zhǔn)確率高，但該方法穩(wěn)定性很差，準(zhǔn)確率波動很大，其海馬體的準(zhǔn)確率相對于平均值波動達(dá)到了10%以上。LW joint方法的準(zhǔn)確率相對于平均值波動在10%以內(nèi)，穩(wěn)定性較好，但只有帶狀體的核團(tuán)分割準(zhǔn)確率超過了MV方法。本文提出的方法提高了三個(gè)核團(tuán)的分割準(zhǔn)確率，尤其是對紋狀體的平均準(zhǔn)確率相對于MV方法提高了8%左右。對帶狀體的分割準(zhǔn)確率提高程度有限，相對于準(zhǔn)確率的平均值波動在5%以內(nèi)，穩(wěn)定性在三個(gè)方法中最高，可以看出本文提出的方法與其他方法相比，有著較好的魯棒性。

表1 三種方法的Dice系數(shù)計(jì)算結(jié)果Table 1 The results of Dice coefficient calculated by the three methods

圖4 Dice系數(shù)Fig.4 Dice coefficient results

圖5為將三個(gè)核團(tuán)的冠狀面分割結(jié)果單獨(dú)呈現(xiàn)的效果圖。第一行至第三行分別為海馬體、紋狀體、屏狀核。第一列至第四列分別為金標(biāo)準(zhǔn)、MV方法分割結(jié)果、LW joint分割結(jié)果、ISRLF分割結(jié)果?？梢钥闯霰痉椒ǖ姆指罱Y(jié)果較其他方法有一定的提升。

圖5 海馬體、紋狀體、屏狀核橫斷面分割結(jié)果圖Fig.5 Transverseplane segmentation results of hippocampus, striatum, and claustrum

4 結(jié)束語

本文提出了一種改進(jìn)的稀疏表示多圖譜分割方法，然后應(yīng)用于獼猴大腦皮下核團(tuán)的分割。該方法以稀疏表示多圖譜分割方法為基礎(chǔ)，將標(biāo)簽信息二值化后引入到圖像塊中，以增加圖像塊的先驗(yàn)信息；然后通過改變信息熵的計(jì)算方式改進(jìn)互信息，將其用來衡量目標(biāo)圖像的目標(biāo)體素所在平面和線性配準(zhǔn)后的圖譜中對應(yīng)平面的全局相似度，使融合時(shí)各圖譜的權(quán)重更加合理。后提出一個(gè)相似性指標(biāo)將非局部塊加權(quán)融合的分割結(jié)果與稀疏表示的融合分割結(jié)果進(jìn)行綜合，從而得到最終的分割結(jié)果。通過與其他方法得到的結(jié)果進(jìn)行Dice系數(shù)的對比，證明本文提出的方法獲得了更高的準(zhǔn)確率和更好的魯棒性。