一種改進稀疏表示的多圖譜腦MR圖像分割算法

曹 敏，侯金秀，張月芳，鄧紅霞，李海芳

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

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

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

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

1 背景描述

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

T=Dα.

(1)

(2)

2 ISRLF方法

2.1 方法框架

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

算法步驟如下：

1) 對圖譜灰度圖像進行線性配準，得到線性配準后的圖像Li.

2) 對線性配準后的圖譜灰度圖像進行非線性配準，得到配準后的結果Fi.將圖譜標簽圖像非線性配準到目標圖像中，得到配準后的結果Ii.

3) 對目標圖像采用局部加權融合的方法進行預分割，得到預分割結果。

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

5) 引入互信息并對信息熵的計算方法進行改進，然后對目標體素和線性配準后的圖譜相對應的平面進行全局的相似性衡量，求平均值，改進稀疏表示計算出的權重。

6) 分別用稀疏表示的標簽融合方法和非局部塊加權的標簽融合方法對目標圖像進行分割，得到各自的分割結果。

7) 采用Dice系數和余弦距離相結合的相似性度量標準對分割結果進行度量，將兩種融合方法的分割結果再次進行融合，最終得到分割結果。

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

2.2 目標圖像標簽塊和圖譜圖像標簽塊的構建

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

本文采用D99-SL、NMT和inia19三個圖譜進行分割。為了將圖譜中的標簽信息引入，需要對目標圖像進行預分割，本方法采用局部加權法對目標圖像進行預分割。構建初始的僅含有灰度信息的圖像塊后，將對應的標簽值塊變成二值化的形式，再與初始圖像塊結合構成新的圖像標簽塊。二值化的具體方式為以中心體素的標簽值為基準，與之標簽值一致的標簽記為1，否則記為0.

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

圖2 圖譜圖像標簽塊的構建Fig.2 Construction of atlas image label block

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

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

2.3 基于互信息的圖像相似度計算

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

(3)

(4)

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

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

(5)

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

(6)

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

(7)

γ=Dice(IAj,IBj) .

(8)

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

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

(9)

(10)

(11)

(12)

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

(13)

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

(14)

W2=W1×S(x) .

(15)

(16)

2.4 標簽融合

基于圖像塊的融合方法有兩種，一種是通過一個相似性度量函數度量所有參與標簽融合的圖像塊與目標圖像塊的相似性的非局部塊加權融合方法，一種是稀疏表示方法，將目標圖像塊由參與融合的圖譜圖像塊稀疏表示。

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

為了結合這兩種方法各自的優(yōu)點，提高對大腦皮下核團的分割準確度，本方法將非局部塊加權方法和稀疏表示的標簽融合方法進行了結合，從而提出了一種新的標簽融合方法ISRLF.

基于圖像塊加權融合方法的過程如下：

1) 將各圖譜灰度圖像和標簽圖像配準到目標圖像T中，得到配準后的灰度圖像Fi和標簽圖像Ii.

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

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

4) 按照公式(18)得到標簽的結果。

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

(17)

(18)

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

(19)

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

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

(20)

3 實驗與結果分析

3.1 數據介紹

實驗選用的測試數據集為牛津大學發(fā)布的獼猴數據集。該數據集用3T的掃描儀采集了20只雄性獼猴的數據。數據集包括T1，靜息態(tài)的fMRI。本實驗使用T1數據。獼猴年齡分布為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，翻轉角度為8°.

3.2 評價指標

實驗選用Dice系數評價分割質量，其計算方法如公式(21)所示。公式中的A表示目標圖像的金標準(即專家分割結果)，B表示用分割方法得出的分割結果?！硎靖骷系捏w素個數，|A∩B|表示集合A與集合B對應體素標簽相同的體素個數，|A|與|B|分別表示集合A與集合B中的體素個數。Dice系數取值范圍為0到1，越接近1，說明準確率越高。

(21)

3.3 結果分析

本實驗選取三個方法分別對海馬體、紋狀體、帶狀體這三個核團的分割結果與金標準進行比較。參與對比的方法有多數投票法(majority voting，MV)，聯合標簽融合法(LW joint)以及本文的方法(ISRLF).

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

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

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

表1 三種方法的Dice系數計算結果Table 1 The results of Dice coefficient calculated by the three methods

圖4 Dice系數Fig.4 Dice coefficient results

圖5為將三個核團的冠狀面分割結果單獨呈現的效果圖。第一行至第三行分別為海馬體、紋狀體、屏狀核。第一列至第四列分別為金標準、MV方法分割結果、LW joint分割結果、ISRLF分割結果?？梢钥闯霰痉椒ǖ姆指罱Y果較其他方法有一定的提升。

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

4 結束語

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