一種新穎的MR脊柱圖像自動分割算法研究

葉偉，陶晶，陳小宇，林千早，魏玲，張小兵，陳萬洪

1.徐州醫科大學附屬淮安醫院 影像科，江蘇 淮安 223002；2.淮安市第四人民醫院 影像科，江蘇 淮安 223002

引言

脊柱是人體中最復雜的承載結構，成人大約長70 cm，分為5個椎體段[1]。椎間盤突出、椎管狹窄和退行性椎間盤是常見的脊柱疾病，MRI是一種行之有效的影像檢查手段，但鑒別病變部位和定量分析往往需要借助圖像分割技術。由于MR圖像自身分辨率和對比度較差，易受噪聲、偽影和局部體積效應等影響，使得MR脊柱圖像自動分割成為一項具有挑戰性的任務。

常用的分割方法有全局閾值法、邊緣檢測法、區域增長方法等[2-4]。灰度閾值由灰度直方圖決定，由于MR圖像的灰度分布具有隨機性，往往導致全局閾值法無法確定最佳分割閾值。此外，閾值分割法忽略圖像像素位置信息，加劇了圖像空間不確定性。邊緣檢測法可用于識別感興趣區域的輪廓邊界，但對邊界的連續性和光滑度要求較高，且需要復雜的后處理操作。區域增長法是閾值法的延伸[5]，將滿足一定灰度相似性測度準則的像素通過鄰域擴展連通成區域，但對初始種子點選取要求較高，且對噪聲敏感。

模糊C均值聚類區別于K均值、C均值等硬性分類算法，是指通過優化目標函數得到每個樣本點對所有類中心的隸屬度，從而決定樣本點的類屬以達到自動對樣本數據進行分類的目的[6]。本文采用基于模糊C均值聚類的分割算法快速準確的提取MR脊柱部分，輔助醫生對脊柱退化狀態分類和畸形的診斷。

1 方法

首先，運用各向異性擴散濾波對原始圖像進行預處理，再采用模糊C均值聚類算法進行圖像分割，最后利用形態學方法進行后處理（圖1）。

圖1 MR脊柱圖像分割流程

1.1 圖像預處理

各向異性濾波是一種基于偏微分方程的濾波技術，建立于熱量的各向異性擴散理論[7]。本研究采用各向異性擴散濾波平滑原始MR脊柱圖像，在圖像平坦區域選擇大尺度平滑，邊緣區域選擇小尺度評卷，能在抑制噪聲的同時保持圖像邊緣信息。各向異性擴散濾波將圖像看作熱量場，每個像素看作熱流，根據當前像素和領域像素關系，判斷是否要向周圍擴散[8]。如果某個鄰域像素和當前像素差別較大，表示此鄰域像素可能處于邊界，則當前像素停止向該方向繼續擴散，確保圖像邊緣完整清晰。主要步驟如下：

首先，設置初始值和邊界值，然后采用公式(1)進行迭代運算得到預處理圖像[9]。其中I表示源圖像，t表示迭代次數。

N、S、E、W代表4個散度，在4個方向上對當前像素求偏導數，由公式(2)給出[10]。

其中N、S、E、W 4個方向的導熱系數由公式(3)給出[11]，其中導熱相關系數k取值越大，平滑效果越好，邊緣保留越差。

根據預實驗結果，綜合考慮計算量和平滑效果，迭代次數取20，導熱系數k取15，λ取為0.15。

1.2 核密度估計的模糊C均值聚類算法優化

模糊C均值聚類算法屬于無監督學習方法，將普通C均值算法對于數據的硬性劃分改進為柔性的模糊劃分[12]。主要思想是使得被劃分到同一簇的對象之間相似度最大，而不同簇之間的相似度最小，已廣泛應用于模式識別和圖像處理研究等領域[13]。模糊C均值聚類算法是通過使得已知數量的聚類目標函數達到最小化來獲得數據集分類，目標函數由公式(4)給出[14]。

其中，m是一個加權指數，0≤m≤∞，通常m取2.0是合理的[15]；uij是群集的隸屬度，0≤uij≤1，且xi是第i個數據點；vj是模糊組j的聚類中心；‖xi-vj‖2是表示數據點與聚類中心間的距離；N表示像素數，而M表示簇中心數，根據預實驗結果取值為3。

模糊聚類是公式(4)所示的目標函數的迭代優化，通過更新隸屬度函數uij和簇中心vj進行。使目標函數達到最小的必要條件由公式(5)給出[16]：

首先選取(0,1)間的隨機數初始化uij，再利用公式(6)計算M個聚類中心vj，和目標函數，如果前后兩次目標函數值的改變量小于某個閾值ε，則迭代停止，否則采用公式(5)更新uij，重復上述操作。ε是迭代停止的臨界值，本文 ε取 0.01。

模糊C均值聚類算法雖然簡潔高效，但初始的M個初始聚類中心vj很難確定，選擇不當，容易陷入局部最優解，本文采用核密度估計來確定模糊C均值聚類所需的M個初始聚類中心，核密度估計是一種非參數估計的建模方法，可直接從圖像的像素連續變化值中估計出概率密度函數，此時圖像像素概率密度估計曲線的峰值點即是初始聚類中心[17]。本文選取高斯函數作為核函數，等概率的選出t個格點的像素灰度值作為觀察值。則高斯核密度估計函數關系式如公式(7)所示[18]：

用公式(8)在圖像中等概率的選出t個格點，窗寬h由公式(9)確定[19]。

其中Ixy為第x行，第y列的像素點灰度值，yj是第j個格點灰度值，max(Ixy)是圖像像素最大值，min(Ixy)是圖像像素灰度最小值，h為窗寬，STD是圖像像素的標準差，IQR為圖像像素的四分位差，N是圖像的像素數量。

1.3 圖像后處理

本研究采用一系列形態學操作從聚類輸出圖像中提取椎體。主要步驟如下：

第一步，先進行2次圖像膨脹，采用3×3的結構元素掃描圖像中的每一個元素，并與其覆蓋的二值圖像進行“與”運算，填充椎體內細小空洞；再進行2次圖像腐蝕，此時結構元素與二值圖像做“或”運算，能夠在不明顯改變物體面積的同時平滑邊界，去除孤立點。

第二步，區域濾波。調用區域濾波函數roifilt2(H,I,BW)對二值圖像BW中灰度值為1（白色）區域進行濾波處理，忽略灰度值為0（黑色）的區域，濾波結果返回到輸出矩陣J中，其中I為目標圖像數據矩陣。

第三步，采用區域測度提取椎體。根據椎體和周圍組織形態分析，椎體縱橫比為1.5～2.0，將此作為區域測度有助于剔除脊柱周圍的韌帶和肌肉。

1.4 圖像質量評價指標

使用 Dice系數（Dice Similarity Coefficient，DC）和HD距離（Hausdorff Distance，HD）評價圖像分割結果，將放射科醫師對圖像的手動分割結果作為參考標準。

Dice系數用于度量兩個集合的一致性[20]，由公式(10)給出。其中X，Y表示兩個集合，分子表示兩個集合的相交操作后的點集合，分母表示兩個集合總和，Dice系數越接近于1，圖像分割結果越精確。

HD距離用于評價兩個二進制點集之間的相似性[21]，見公式(11)。其中，h(X,Y)稱為從集合X～Y的單向HD，見公式(12)。HD是測量點的相似位置，能夠抵抗點位置的微小波動，準確測量點集之間非相似性。HD越接近于0，圖像分割越精確。

2 結果與討論

為了定性、定量評估本文提出算法的優越性，強健性和實用性，選取5幅MR脊柱T1加權矢狀位平面圖像作為實驗對象，圖像矩陣大小均為512×512，具有512個灰階，所用算法均在MATLAB 2013a編程環境下仿真實現。

2.1 不同算法脊柱分割結果

應基于模糊C均值分割算法、Otsu閾值分割算法和K均值分割算法的MR脊柱提取過程，見圖2～4。首先對源圖像進行預處理，接著采用不同算法進行分割，然后采用腐蝕、濾波和基于縱橫比的形態學處理。結果表明，模糊C均值聚類分割方法能完整的保留脊柱信息，而基于Otsu閾值算法和K均值聚類方法均只分割出部分椎體，脊柱信息嚴重不完整，這是由于Otsu法難以確定最佳閾值，導致分割不足；K均值算法分聚類中心和數目對分割結果影響很大，容易出現分割不足和過分割。

2.2 臨床實例分析

為了驗證本文分割算法的臨床應用價值，選取4例輕度腰痛患者的矢狀位MR脊柱T1加權圖像，其中包括2名女性和2名男性，年齡為45～60歲。4名患者的脊椎分割結果，見圖5，由模糊C均值聚類分割后進行標記著色與原始圖像疊加所得。可見模糊C均值聚類分割算法能夠適用于多種不同的脊柱類型，為臨床醫生診斷提供便利。

圖2 基于模糊C均值聚類分割結果

圖3 基于Otsu閾值分割結果

圖4 基于K-均值聚類分割結果

圖5 4例腰痛患者脊柱分割結果

此外，分別采用Otsu閾值算法、K均值聚類算法和模糊C均值聚類算法對4例患者的脊柱圖像進行分割，并以放射科主任醫師手動分割結果作為驗證參考標準，分別計算不同分割算法所得的DC值和HD值，結果見圖6。DC值大小排序依次為模糊C均值聚類算法、K均值分割算法和Otsu閾值法，其中基于本文提出算法的DC值分別為0.835，0.900，0.852和0.830（接近于1），HD值在模糊C均值分割算法中最小，分別為0.040，0.040，0.036和0.040，其次是K均值聚類分割算法，Otsu閾值分割算法的HD值最大。表明基于模糊C均值分割算法所得的分割結果與參考標準相似性最強，分割結果最穩定，其次是K均值聚類算法，Otsu閾值分割效果最差。

圖6 不同分割算法的Dice系數（a）和HD距離（b）比較

3 結論

本文提出了一種基于改進的模糊C均值聚類的無監督自動分割算法，并應用于脊柱MR冠狀位平面圖像的分割。其中各向異性濾波函數能平滑原始圖像，高斯核密度估計可以為模糊C均值聚類算法選出最佳的聚類的中心和聚類數目，模糊C均值算法得到二值化初始分割圖像，一系列的形態學操作成功提取出椎體感興趣區域，對于退行性腰椎病變等腰椎疾病的定量分析具有重要指導意義。定量和定性評估均表明模糊C均值聚類算法能快速準確提取出脊椎區域，明顯優于K均值聚類算法和Otsu閾值法，具有很強的魯棒性，是一種可行的MR脊柱圖像分割算法。

[參考文獻]

[1] 劉松濤,殷福亮.基于圖割的圖像分割方法及其新進展[J].自動化學報,2012,38(6):911-922.

[2] 王江濤,練煜,石紅巖,等.結合模糊C均值聚類和邊緣檢測算法的彩色圖像分割[J].蘭州文理學院學報(自然科學版),2015,29(2):61-65.

[3] 劉俊,馬燕,陳坤.一種改進的基于區域生長的彩色圖像分割方法[J].計算機應用與軟件,2012(11):288-291.

[4] 曾文權,何擁軍,崔曉坤.基于各向異性濾波和空間FCM的MRI圖像分割方法[J].計算機應用研究,2014,31(1):316-320.

[5] 付麗娟,姚宇,付忠良.中值濾波與各向異性擴散相結合的醫學圖像濾波方法[J].計算機應用,2014,34(1):145-148.

[6] Kannan SR,Ramathilagam S,Sathya A,et al.Effective fuzzy c-means based kernel function in segmenting medical images[J].Comput Biol Med,2010,40(6):572-579.

[7] Li BN,Chui CK,Chang S,et al.Integrating spatial fuzzy clustering with level set methods for automated medical image segmentation[J].Comput Biol Med,2011,41(1):1-10.

[8] Gong M,Liang Y,Shi J,et al.Fuzzy C-means clustering with local information and kernel metric for image segmentation[A].IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society[C].2013,22(2):573-584.

[9] 楊素華,陳瓊,羅艷芬.基于Graph-Cuts的腦部MRI圖像腦組織提取方法[J].中國生物醫學工程學報,2014,33(5):525-531.

[10] 黃會營,莊淑君.改進的正交匹配追蹤圖像快速重建算法[J].光學技術,2014,40(6):515-519.

[11] 謝永華,徐趙飛,范文曉.基于高斯尺度空間粗糙度描述子的花粉圖像分類識別[J].計算機應用,2015,35(7):2039-2042.

[12] 陳金位,吳冰.二維直方圖重建和降維的Otsu閾值分割算法[J].圖學學報,2015,36(4):570-575.

[13] 徐超,黃風華,毛政元.一種改進的二維Otsu閾值分割算法[J].電子技術應用,2016,42(12):108-111.

[14] Ananthi VP,Balasubramaniam P,Lim CP.Segmentation of gray scale image based on intuitionistic fuzzy sets constructed from several membership functions[J].Pattern Recogn,2014,47(12):3870-3880.

[15] 周嘯虎,高偉,張子齊.基于質子密度和弛豫時間的大腦MR圖像分割新算法[J].中國醫療設備,2016,31(10):25-28.

[16] 宋國權,李金鋒.基于聚類算法的腦部MR圖像分割[J].中國醫療設備,2017,32(1):26-29.

[17] Gómez D,Yá?ez J,Guada C,et al.Fuzzy image segmentation based upon hierarchical clustering[J].Knowl-Based Syst,2015,87(C):26-37.

[18] Dubey YK,Mushrif MM.FCM clustering algorithms for segmentation of brain MR images[J].Adv Fuzzy Syst,2016(3):1-14.

[19] 朱然,李積英.幾種優化FCM算法聚類中心的方法對比及仿真[J].計算機技術與發展,2015(5):17-20.

[20] ?ajn L,Kononenko I,Mil?inski M.Computerized segmentation and diagnostics of whole-body bone scintigrams[J].Comput Med Imag Grap,2007,31(7):531-541.

[21] 李俊林,符紅光.改進的基于核密度估計的數據分類算法[J].控制與決策,2010,25(4):507-514.