基于模糊隸屬度函數的腦圖像分割方法

楊 兵 劉曉芳

(中國計量大學計算機應用與技術研究所 浙江 杭州 310018) (浙江省電磁波信息技術與計量檢測重點實驗室 浙江 杭州 310018)

0 引 言

磁共振成像(Magnetic Resonance Imaging,MRI)具有軟組織對比度高、無放射損傷等優點[1]，廣泛用于人體腦部成像。腦圖像分析是了解人腦健康狀況以及生理特征的重要渠道，其中腦圖像分割是研究人腦解剖結構以及智能診斷腦部疾病的基礎，近年來，機器學習技術因其具有良好的數據預測能力以及模型泛化性能，被廣泛用于醫療數據智能診斷[2]中。腦組織分割和腦腫瘤智能檢測是腦部圖像分割的兩個重要研究方向，本文的目的是應用機器學習技術對腦圖像中腦灰質(Gray Matter,GM)、腦白質(White Matter,WM)、腦脊液(Cerebrospinal Fluid,CSF)三種主要腦組織進行自動分割。

由于成像設備等原因，磁共振圖像存在各組織之間沒有明顯的邊界、各組織之間灰度值重疊等現象[1]，針對此問題，Vaishnavee等[3]將支持向量機用于腦影像智能分析中，證明支持向量機可以根據醫學圖像灰度不均勻以及組織邊緣不清晰等特點做出準確的識別分析。支持向量機是一種基于類間最大距離的數據處理算法，其分類超平面的建立僅依賴于少數支持向量[4]，具有較好的泛化性能和判別能力，在此基礎上，模糊支持向量機將模糊理論和支持向量機相結合[5]，提出根據空間距離判別規則對樣本賦權，在解決樣本數據不完全屬于某一類的情況下表現出了較好的預測能力，廣泛應用于醫學圖像分析中。

隸屬度函數的設計是構建模糊支持向量機模型的關鍵[5]。目前，隸屬度函數的設計大多基于同類樣本之間的幾何距離，沒有考慮到不同類樣本之間的模糊距離，無法有效解決灰度不均勻帶來的像素值交叉重疊問題?；诖?，本文提出一種基于類間樣本模糊距離的隸屬度確定方法，綜合考慮同類樣本中心以及不同類樣本中心對待賦權樣本的影響，訓練模糊支持向量機模型對三種主要腦組織進行自動分割。

1 相關工作

1.1 圖像分割

圖像分割是指利用計算機輔助技術將圖像劃分為具有不同性質的圖像塊[6]，傳統算法如直方圖方法[7]提取圖像特征如梯度特征，顏色直方圖、紋理特征將圖像劃分為不同圖像塊。然而，由于醫學圖像具有灰度不均勻、各組織之間邊界不清晰等特點，傳統算法不能較好適應此種情況。

研究人員將支持向量機應用于醫學圖像處理中，如Blumenthal等[8]將支持向量機應用于腫瘤圖像分類，Khedher等[9]將支持向量機與主成分分析方法結合用于阿爾茲海默癥的早期診斷，Emblem等[10]將支持向量機用于術前腦部膠質瘤生存預測。然而，上述方法沒有充分考慮到醫學圖像的特點，沒有將同類樣本距離以及不同類樣本距離引入支持向量機中，無法有效解決醫學圖像分割面臨的邊界不清晰和灰度不均勻等問題。由支持向量機改進的模糊支持向量機在處理樣本不完全屬于某一類的問題上表現出不錯的應用前景，模糊支持向量機的構建主要分為三部分：1) 根據樣本隸屬度函數計算樣本隸屬度；2) 對樣本進行賦權，訓練模糊支持向量機模型；3) 對模型調整參數，使其具有較好的判別能力與泛化性能。

1.2 支持向量機

支持向量機是一種基于類間距離最大化的數據處理算法[11]，其目標是在多維特征空間中確定最優分類超平面，使得分類確信度(confidence)最大。支持向量機能夠很好地處理線性可分問題，對于線性不可分問題，支持向量機引入核函數將原始數據從低維空間映射到高維空間，從而解決原始數據在低維空間中線性不可分的問題。即考慮以下問題：

(1)

式中:φ為原始數據空間到高維特征空間的映射；w為數據權重項；b為數據偏置項。核函數就是將映射形式φ隱式地表示為原始數據內積的形式：

(2)

(3)

s.t.αi≥0,i=1,2,…,n

對于離群點，支持向量機通過引入松弛變量ξ來更好地容忍少數離群點對于分類超平面地影響，目標是解決以下原始問題：

(4)

s.t.yl(wTxl+b)≥1-ξl,i=1,2,…,n

ξl≥0,l=1,2,…,n

式中：w為權重因子；C用于控制最大分離超平面與數據偏差的權重；x為樣本數據點；ξ為松弛變量；y為樣本數據的標簽。通過拉格朗日法求解得到以下對偶問題：

(5)

0≤αi≤C，i=1,2,…,n

通過對拉格朗日乘子α的求解，將原始問題轉化為對偶問題。

1.3 模糊支持向量機

傳統支持向量機中，每個樣本點完全屬于某一類，不能處理樣本同時屬于多類的情況。Khodagholi等[12]將模糊理論引入支持向量機中，賦予每個樣本點不同的隸屬度ml(0

s.t.yl(w·xl+b)≥1-ξl

ξl≥0,l=1,2,…,n

(6)

由拉格朗日乘子法得到原始問題的對偶問題：

0≤αj≤mjC,j=1,2,…,n

(7)

樣本隸屬度的確定是構建模糊支持向量機的關鍵。Lin等[13]提出一種基于樣本到樣本中心距離的隸屬度確定方法，假設有樣本集{(x1,y1),(x2,y2),…,(xn,yn)}，y∈{-1,+1}，x+表示正類樣本的均值，x-表示負類樣本的均值，r+=max|xi-x+|，r-=max|xi-x-|，則樣本隸屬度mi定義為：

(8)

上述方法沒有考慮到各類別樣本之間的內在聯系，將各類別的樣本孤立看待。張翔等[14]考慮了樣本與樣本中心的關系以及樣本之間的緊密度，提出了一種基于樣本緊密度的隸屬度計算方法，但無法有效解決孤立樣本隸屬度賦權問題。

2 模糊隸屬度函數設計

由于成像技術等原因，醫學圖像中的每個像素的灰度值通常是多種組織灰度值的平均，這也導致醫學圖像中組織與組織之間邊界不清晰。針對此問題，本文設計了一種基于不同類別樣本之間模糊距離的隸屬度計算方法。

圖1中，超球體分別表示三種不同的腦組織。AB、AC、BC表示各腦組織樣本間最遠距離。a、b、c分別表示各腦組織樣本中心，設樣本點x1到三種腦組織的距離分別為L2、L1、L3。

圖1 樣本隸屬度

(9)

樣本點x2到超球體中心b、a、c的距離分別為M1、M2、M3。

(10)

樣本點x3到超球體中心b、a、c的距離分別為N1、N2、N3。

(11)

則樣本點x1的隸屬度m1為：

(12)

樣本點x2的隸屬度m2為：

(13)

由于樣本點x3位于三角形ABC之外，故樣本點x3的隸屬度m3為：

(14)

式中：λ、μ為模糊距離系數。若L1=M1，即樣本點x1和樣本點x2在超球體內處于同一球面上，但由于L2+L3m2，這表明，即使處于同一超球面上的樣本點，距離其他兩類中心距離更近的樣本點獲得更大的隸屬度。若有L1=N1，即樣本點x1和樣本點x3在超球體內處于同一球面上，但由于樣本點x3位于三角形ABC之外，故樣本點x3賦予較小的隸屬度，從而使得位于三類樣本中心附近的樣本點獲得更大隸屬度。

綜上，對于位于三角形ABC之內的樣本隸屬度計算公式為：

(15)

(16)

式中：μ為模糊距離系數，使位于三角形ABC之外的樣本隸屬度mi<0.5。

3 整體算法流程

本文提出的算法整體流程如圖2所示。

根據圖2，模糊隸屬度函數算法主要步驟如下：

1) 根據具體應用實例，提取數據特征，如圖像的梯度特征，紋理特征等。

2) 結合樣本標簽，計算每個數據樣本的模糊隸屬度。應用于圖像數據時，不計算背景樣本的模糊隸屬度。

3) 根據輸入數據特征，并劃分數據集為訓練集、驗證集以及測試集，將步驟2)計算的樣本模糊隸屬度作為樣本數據權重，建立模糊支持向量機模型。

4) 使用交叉驗證策略，根據驗證集的分類精度調整模型參數(懲罰系數C、核函數系數、模糊距離系數μ和λ等)。

5) 使用步驟4)得到的分類模型，驗證模型在測試集上的精度。

此外，本文還給出了具體的算法執行流程，其中圖像的問題特征提取和樣本的模糊隸屬度計算以偽代碼方式展現。

Step1設定模型初始參數，包括懲罰系數C，迭代次數MaxG，核函數系數σ，模糊距離系數λ、μ等。

Step2數據樣本預處理，包括歸一化，去除異常值等。

(17)

Step3設定數據劃分比例：a、b、c。分別將數據分別劃分為訓練集、驗證集、測試集。

Step4提取數據特征，如圖像的梯度特征、紋理特征(對比度CON、ASM能量、熵)等。

Gx(x,y)=H(x+1,y)-H(x-1,y)

Gy(x,y)=H(x,y+1)-H(x,y-1)

(18)

以圖像的紋理特征提取為例，imgLength、imgwidth分別表示圖像的長和寬，GLCM_GrayLevel表示灰度共生矩陣的灰度級數，greycomatrix_diatance表示灰度共生矩陣距離，greycomatrix_angles表示灰度共生矩陣的角度，texture_property表示需要提取的紋理特征，glcm表示灰度共生矩陣。

For i=1 to imgLength:

For j=1 to imgWidth:

If i<3 or j<3:

continue

# 使用7×7的窗口大小計算灰度共生矩陣

GLCM_window=img[i-3:i+4,i-3,j+4]

# 計算灰度共生矩陣

glcm=

greycomatrix(glcm_window,

greycomatrix_diatance,greycomatrix_angles,GLCM_GrayLevel)

# 計算紋理特征1

texture_property_1=

greycoprops(glcm,texture_property[0])

texture_1[i,j]=texture_property_1[0][0]

# 計算紋理特征2

texture_property_2=

greycoprops(glcm,texture_property[1])

texture_2[i,j] = texture_property_2[0][0]

# 計算紋理特征3

texture_property_3=

greycoprops(glcm, texture_property[2])

texture_3[i,j]=texture_property_3[0][0]

Step5計算并保存每個數據樣本的模糊隸屬度：

(19)

其中，兩類樣本之間的最遠距離用動態規劃法求解，N表示負類樣本數，M表示正類樣本數，Dn(i,j)表示負樣本i到正樣本j的距離，Dp(i,j)表示正樣本i到負樣本j的距離，偽代碼如下：

For i=1 to N:

FarthestDistance[i]=0

For j=1 to M:

FarthestDistance[j]=0

For i=1 to N:

For j=1 to M:

If Dn(i,j)+Dp(i,j)>FarthestDistance[i][j]:

FarthestDistance[i][j]=Dn(i,j)+Dp(i,j)

If Dn(i,j)+Dp(j,i)>FarthestDistance[i][j]:

FarthestDistance[i][j]=Dn(i,j)+Dp(j,i)

Step6選擇并初始化核函數、懲罰系數C、模糊距離系數λ和μ等。

Step7根據Step 1-Step 6建立模糊支持向量機模型。

(20)

0≤αj≤mjC,j=1,2,…,n

Step8選擇評價指標對驗證集定律評價。

Step9若Step 8的評價結果不滿足精度要求或算法迭代數小于迭代次數MaxG，則返回Step 6，否則執行Step 10。

Step10保存算法模型，輸入測試集并輸出預測結果，算法結束。

4 實 驗

實驗采用的硬件設備是Intel Core 4.2 GHz i7- 7700k CPU以及NVIDIA Geforce GTX 1080Ti GPU，軟件環境采用Python 3.5，所用圖像處理軟件包括scikit-image，SimpleITK、以及神經影響處理工具NiBabel。此外，還用到機器學習處理軟件包scikit-learn(sklearn)。

4.1 實驗數據

本文所用實驗數據來自brainweb[15]，其包含了多個數據的人工分割結果。磁共振圖像中縱向弛豫時間T1是指成像掃描儀射頻脈沖發射后，核自旋磁化強度受到激發時，從高能狀態恢復到原始穩定狀態的快慢程度[16]。T1值也是磁共振成像中軟組織對比的重要來源，針對不同組織的成像，T1映射圖可以作為反映組織生理健康程度的參考，故計算T1映射圖對于臨床診斷、疾病治療都至關重要。

本文根據T1值生成圖像標簽。在成像序列中，重復時間TR、回波時間TE、反轉時間(Inversion Time，TI)以及發射角(Flip Angle，FA)影響圖像對比度。根據成像對比度質量和實際操作經驗設置SE序列以及FLASH序列的成像參數。此外，層面厚度ST由梯度強度和脈沖帶寬決定，在2D圖像中，ST越小則空間分辨率越高，但圖像信噪比也隨之下降，本文結合磁共振相關原理和常用層面厚度來設置層面厚度ST。磁共振圖像序列采用自旋回波-反轉恢復(Spin-Echo Inversion Recovery,SE-IR)序列，序列翻轉時間(Inversion Time,TI)分別為800、1 200、1 600、2 200采集得到，如圖3所示。圖3成像參數為：重復時間(Repetition Time,TR)/回波時間(Echo Time,TE)為2 500/14，層面厚度(Slice Thickness,ST)為3 mm。

(a) TI=800(b) TI=1 200 (c) TI=1 600 (d) TI=2 200圖3 不同TI的T1加權圖像序列

待分割腦圖像如圖4所示，實驗中采用了兩個脈沖序列：SE序列和小角度快速激發梯度回波(Fast Imaging using Low Angle Shot,FLASH)序列。圖像的噪聲水平為3%和20%，灰度不均勻性(Intensity Non-uniformity,INU)為20%和40%。其中SE成像參數：TR/TE為1 500/13，ST為3 mm。FLASH成像參數：TR/TE為25/13，發射角(Flip Angle)為15°。

(a) SE圖像，從左至右噪聲比例和灰度不均勻程度 依次為：3%、20%，3%、40%，20%，20%

(a) FLASH圖像，從左至右噪聲比例和灰度不均勻程度 依次為：3%、20%，3%、40%，20%，20%圖4 部分實驗圖像

4.2 標簽圖像獲取

本文將T1映射圖作為磁共振圖像數據的參考標簽。首先計算MR掃描部位的T1值，通過設置脈沖序列中特定的成像參數，掃描得到不同參數MR圖像序列，再根據Bloch方程設計T1的計算方法[17]。最后，通過非線性最小二乘法，擬合得到每個像素點的T1值。本文中，采用經典的自旋回波反轉恢復序列，通過采集不同反轉時間的T1加權圖像，然后運用計算T1值：

(21)

式中：a、b是待定復數形式參數，通過擬合得到；TI是反轉時間(Inversion Time)，是采用自旋回波反轉恢復序列進行掃描時需要設置的參數；n=1,2,…是序列圖像的序號；S表示采集的MR圖像上的某個像素值。

4.3 標簽圖像處理

為了從T1映射圖上得到訓練分割模型所需要的標簽信息，需要對T1映射圖做一些預處理。圖5為原始T1映射圖的直方圖，可以看出，每種腦組織像素值在直方圖上表現為一個數值范圍(其中0表示背景)，對其取三個閾值再經過形態學濾波去除圖像空洞及毛刺得到模糊標簽。由于模糊標簽中存在一些噪聲標簽，影響后續腦組織分割，在本文中，采用Northcutt等[18]提出的標簽去噪方法去除噪聲標簽，最終得到的模糊標簽為分割腦組織的圖像標簽。

圖5 原始T1映射圖的灰度直方圖 圖6 原始T1 映射圖 圖7 經過預處理后的T1映射圖

4.4 模型評價指標

采用的分割評價指標包括分類準確率(Classification Accuracy，CA)、召回率(Recall)、Dice系數(Dice Index，DI)。三者計算公式如下：

(22)

(23)

(24)

式中：a表示正確分類樣本數；b為錯誤分類樣本數；TP表示分類器將正類劃分為正類的樣本數；FP表示分類器將負類劃分為正類的樣本數；FN表示分類器將正類劃分為負類的樣本數。

4.5 腦圖像分割實驗

本文評估了不同組合的圖像特征以及模型參數(模糊支持向量機懲罰系數C、高斯核函數寬度、模糊距離系數μ和λ)對SE圖像中三種腦組織分割結果的影響，采用的圖像為3%noise、20% INU。模糊距離參數λ和μ對分類精度的影響見表1，懲罰系數C和核函數寬度σ對分類精度的影響見表2，不同特征組合對分類精度的影響見表3。

表1 模糊距離參數λ和μ對分類精度的影響

表2 懲罰系數C和核函數寬度σ對分類精度的影響

續表2

表3 不同特征組合對分類精度的影響

特征組合1：灰度共生矩陣窗口大小為5×5，像素距離1 mm，像素角度45度，梯度，對比度(contrast)，相關性(correlation)，ASM能量(Angular Second Moment,ASM)，逆差矩(Inverse Different Moment,IDM)。

特征組合2：灰度共生矩陣窗口大小為5×5，像素距離3 mm，像素角度90度，梯度，相關性(correlation)，ASM能量(Angular Second Moment,ASM)。

特征組合3：灰度共生矩陣窗口大小為5×5，像素距離3 mm，像素角度45度，梯度，逆差矩(Inverse Different Moment,IDM)，ASM能量(Angular Second Moment,ASM)。

特征組合4：灰度共生矩陣窗口大小為7×7，像素距離3 mm，像素角度90度，梯度，對比度(contrast)，相關性(correlation)，ASM能量(Angular Second Moment,ASM)。

特征組合5：灰度共生矩陣窗口大小為7×7，像素距離3 mm，像素角度90度，像素角度45度，梯度，對比度(contrast)，相關性(correlation)，ASM能量(Angular Second Moment,ASM)，逆差矩(Inverse Different Moment,IDM)。

特征組合6：灰度共生矩陣窗口大小為7×7，像素距離3 mm，像素角度90度，像素角度45度，對比度(contrast)，相關性(correlation)，ASM能量(Angular Second Moment,ASM)，逆差矩(Inverse Different Moment,IDM)。

本文基于灰度共生矩陣(Gray Level Co-occurrence Matrix,GLCM)提取腦組織特征。GLCM是一種在特定空間和方向提取像素間相關性的紋理描述方法，由于圖像像素值分布導致圖像空間像素出現一定的排列規律，故圖像中間隔一定距離的像素之間存在可提取的紋理特征?；诒?-表3給出的不同模型參數對分類精度的影響，并結合腦部圖像呈現細小的紋理性質，本文基于GLCM，用7×7的滑動窗口，像素距離1 mm和3 mm，像素角度45度和90度提取了腦圖像的紋理特征，其統計量分別為：對比度(contrast)，相關性(correlation)，逆差矩(Inverse Different Moment,IDM)，ASM能量(Angular Second Moment,ASM)，共提取紋理特征總計16個，此外，還提取了梯度值等特征10個。

本文所用腦圖像均重采樣為60×180×180大小，為了增加圖像數量，還使用了旋轉、中心裁剪等數據增強技術，得到擴充后圖像8 000幅，其中，5 000幅用于模型訓練，1 000幅用于模型驗證，2 000幅用于模型測試。模糊支持向量機中，懲罰系數C=1 000，高斯核函數寬度σ=1,模糊距離系數λ、μ分別取0.6、0.4。本文分別在SE圖像和FLASH圖像，不同噪聲比例和灰度不均勻程度上進行了實驗，分割結果如圖8所示。

(a) 標準圖像(Ground Truth)及其標準分割結果

(b) SE圖像及腦組織分割結果(腦脊液，腦灰質，腦白質)

(c) FLASH圖像及腦組織分割結果(腦脊液，腦灰質，腦白質)圖8 各腦組織分割結果

可以看出，在SE圖像上，噪聲3%，灰度不均勻40%的圖像的腦白質分割效果最好，也最接近標準圖像的分割結果。本文方法能較好地判別孤立像素點，分割結果空洞情況較少，能對不同序列圖像做出準確分割。

從表4中可以看出，對于SE圖像，以CA為評價指標時，本文方法對噪聲3%的圖像實現了精確分割。結合表4和表5，SE圖像的分割結果稍好于FLASH圖像，在噪聲20%、灰度不均勻20%圖像上，兩者相差較多。在SE圖像上，本文方法平均CA為87%，平均召回率為0.81，平均DI系數為0.82。

表4 SE圖像分割評價

表5 FLASH圖像分割評價

上述實驗結果表明，本文方法對于不同序列的腦圖像具有較好的分割準確率和泛化性能。此外，本文還給出了三種腦組織在三種評價指標下的分割精度，如表6-表8所示。

表6 腦脊液分割評價

表7 腦白質分割評價

表8 腦灰質分割評價

可以看出，腦白質的分割精度最高，對腦白質的分割誤判也最少，其次是腦脊液，最后是腦灰質。三組圖像中腦白質的平均分割準確率達到90.33%，腦脊液平均分割準確率為83%，腦灰質平均分割準確率為83%。在實驗樣本圖像中，腦白質的組織邊緣較清晰，基于灰度共生矩陣提取腦組織的紋理特性時，腦白質的紋理可以較好地提取出來。

此外，為了說明本文分割方法在腦組織分割上的有效性與準確性，本文方法與目前主流的腦組織分割方法隨機森林(Random Forest)方法、模糊C均值方法(Fuzzy C means)、Snake形變模型進行分析比較，結果如表9和圖9所示。

表9 本文方法與其他三種方法分割性能比較

表9中，對于噪聲3%、灰度不均勻20%的圖像，本文方法DI系數為0.90，高于其他三種方法中最高的Snake形變模型0.89的DI系數，相對于隨機森林分割算法0.85的DI系數高出5.9%。對于低噪聲比例和低灰度不均勻程度的圖像來說，本文方法和模糊C均值所取得的分割精度接近。對于灰度不均勻為40%的圖像，本文方法DI系數為0.87，而隨機森林算法所取得的DI系數為0.83。對于噪聲20%的圖像，相比于其他算法，本文方法的DI系數明顯提升，和隨機森林比較，本文方法DI提升10.4%，相對于Snake形變模型和模糊C均值分別提升4.2%、7.2%。圖9表明，本文方法相對于其他三種方法在兩種不同的Noise比例和INU的圖像上的分割精度均有所提升。

5 結 語

針對醫學圖像組織邊緣不清晰以及灰度不均勻問題，本文提出了一種基于類間模糊距離的隸屬度函數計算方法，樣本隸屬度的確定不僅取決于同類樣本之間的空間距離，還取決于不同類樣本之間的空間距離，以解決樣本標簽交叉重疊的問題。此外，針對醫學圖像人工標注的困難，本文提出了基于T1映射圖的模糊標簽獲取方法，減少人工標注圖像的困難。三種主要腦組織分割實驗表明，本文方法可準確有效地分割腦部圖像，在處理醫學圖像具有很好應用前景。此外，本文提出的模糊隸屬度函數在模糊支持向量機的基礎上，可用于分類噪聲較多的數據樣本。一方面，模糊支持向量機對通過樣本重新調整權重的方式增加了對噪聲樣本的魯棒性；另一方面，基于本文提出的模糊隸屬度函數在綜合考慮類內距離和類間距離的基礎上，重新計算樣本的模糊隸屬度，進一步提升了算法模型對噪聲樣本的魯棒性。