基于Nystrom方法的水平集醫(yī)學(xué)圖像分割算法

王曉丹，張龍波，王 雷，劉 晨

(山東理工大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，山東 淄博 255049)

1 相關(guān)研究

近年來(lái)，水平集算法已經(jīng)較為成熟地應(yīng)用于醫(yī)學(xué)圖像分割領(lǐng)域，這是一種有效解決曲線演化問(wèn)題的數(shù)值方法，并且計(jì)算穩(wěn)定，適宜任意位數(shù)空間[1].水平集的基本思想是低維到高維的映射，把低維空間上的函數(shù)通過(guò)水平集的方法轉(zhuǎn)化到高維空間.文獻(xiàn)[2]提出了著名的測(cè)地線主動(dòng)輪廓模型(Geodesic Active Contour，GAC)，它允許曲線在演化過(guò)程中改變拓?fù)浣Y(jié)構(gòu)，分割結(jié)果與初始化曲線的拓?fù)湫螒B(tài)無(wú)關(guān).但需要對(duì)水平集函數(shù)不斷初始化，代價(jià)昂貴，且演化曲線易在弱邊緣處泄露，造成分割錯(cuò)誤.

由于圖像分割問(wèn)題是將圖像的像素進(jìn)行分類，因此可以運(yùn)用聚類分析優(yōu)化圖像分割算法.文獻(xiàn)[3]提出基于模糊聚類水平集的圖像分割方法，算法首先使用模糊C均值(fuzzy C-means clustering，F(xiàn)CM)聚類，得到聚類后的圖像模型，再利用水平集算法進(jìn)行圖像分割.雖然提高了分割的精度，有效地減少了迭代次數(shù)，但FCM方法對(duì)初始化參數(shù)敏感，容易陷入局部最小值.文獻(xiàn)[4]提出了基于核模糊聚類的變分水平集醫(yī)學(xué)圖像分割算法，將原始圖像進(jìn)行核模糊C均值聚類，然后將聚類結(jié)果帶入水平集函數(shù)得到初始輪廓，最后利用一種無(wú)須重新初始化的水平集模型分割圖像.核函數(shù)對(duì)噪聲有很強(qiáng)的處理能力，但基于核方法的模糊聚類算法迭代次數(shù)多，運(yùn)算時(shí)間長(zhǎng).除了模糊聚類，基于圖論的聚類也是近年來(lái)研究的熱點(diǎn).文獻(xiàn)[5]介紹并使用基于連續(xù)最大流的方法對(duì)心臟CT圖像進(jìn)行分割實(shí)驗(yàn)，最終得到左心室心機(jī)的分割結(jié)果.文獻(xiàn)[6]將圖論方法應(yīng)用于脊椎骨磁共振圖像的分割.文獻(xiàn)[7]提出一種新的基于圖論的分割方法，將圖像中灰度相同的像素作為一類，改進(jìn)加權(quán)函數(shù)的定義，將節(jié)點(diǎn)與區(qū)域間的空間近鄰關(guān)系約束進(jìn)加權(quán)函數(shù).

針對(duì)現(xiàn)有的基于最小化區(qū)域擴(kuò)展擬合能量的分割模型對(duì)邊緣模糊、噪聲強(qiáng)的圖像迭代時(shí)間久，易產(chǎn)生邊緣泄露的問(wèn)題，本文提出一種基于Nystrom方法的水平集醫(yī)學(xué)圖像分割算法.

2 基于圖論的圖像分割

圖論法要求將圖像映射為一個(gè)帶權(quán)的無(wú)向圖，把像素當(dāng)作節(jié)點(diǎn)，構(gòu)成最小支撐樹(shù)，利用最小剪切準(zhǔn)則得到閾值從而實(shí)現(xiàn)圖像的最優(yōu)分割[8].

2.1 譜分割

譜聚類算法是將聚類問(wèn)題轉(zhuǎn)化為一個(gè)無(wú)向帶權(quán)圖的多路劃分問(wèn)題[9].無(wú)向圖的頂點(diǎn)V={v1,v2,…vn}表示像素集，vi和vj之間的相似度權(quán)值Wij表示無(wú)向圖的邊E，將聚類問(wèn)題轉(zhuǎn)化為尋找最優(yōu)的分割準(zhǔn)則G(V,E)劃分為子圖內(nèi)部相似度最大，子圖之間相似度最小.經(jīng)典的NCut(Normalize cut)方法常作為最優(yōu)分割準(zhǔn)則，定義圖G劃分的兩個(gè)子圖A∪B=V，A∩B=?，代價(jià)函數(shù)為

(1)

為了實(shí)現(xiàn)最佳二分割，必須找到一個(gè)合適的A和B集合使代價(jià)函數(shù)最小.根據(jù)圖譜理論，可以通過(guò)計(jì)算拉普拉斯矩陣的第二最小特征值對(duì)應(yīng)的特征向量對(duì)圖進(jìn)行分割.這樣就把圖劃分的問(wèn)題轉(zhuǎn)化為求解矩陣特征值和特征向量[10].本文使用高斯核函數(shù)度量像素之間的相似度

2.2 基于Nystrom方法的譜分割

譜分割在理論和應(yīng)用方面效果都不錯(cuò)，但仍有許多不足之處.隨著圖像的分辨率增大，相似矩陣成冪級(jí)數(shù)增長(zhǎng)，花費(fèi)大量的存儲(chǔ)空間和計(jì)算時(shí)間，F(xiàn)owlkes[11].提出通過(guò)Nystrom采樣可近似估算相似矩陣和特征向量. 設(shè)原始圖像有N個(gè)像素，從中隨機(jī)選取m個(gè)像素點(diǎn)，剩下n=N-m個(gè)像素點(diǎn)m?N，構(gòu)造相似矩陣W為

(2)

3 本文算法

本文提出的基于Nystrom方法的水平集醫(yī)學(xué)圖像分割算法，首先利用Nystrom譜聚類的思想，對(duì)原始圖像隨機(jī)采樣，然后分別計(jì)算樣本點(diǎn)之間的距離和樣本點(diǎn)與非樣本點(diǎn)之間的距離，并將距離轉(zhuǎn)化為相似矩陣，拉普拉斯歸一化處理后，得到k個(gè)特征值及對(duì)應(yīng)的特征向量，再利用k-means聚類特征向量，得到聚類的圖像，最后利用Li提出的水平集分割方法[13]分割圖像得到結(jié)果.具體步驟為：

(1)輸入原始圖像，通過(guò)Nystrom算法對(duì)像素點(diǎn)隨機(jī)采樣，選取一定數(shù)量的樣本點(diǎn).計(jì)算樣本點(diǎn)之間的歐式距離、樣本點(diǎn)和非樣本點(diǎn)之間的歐氏距離.

(2)將距離轉(zhuǎn)化為相似矩陣A和B，拉普拉斯歸一化相似矩陣.執(zhí)行正交化和特征分解，得到k個(gè)特征值和對(duì)應(yīng)的特征向量，并按降序排列.

(3)利用k-means算法聚類，得到聚類后的圖像.

(4)初始化水平集函數(shù)、高斯核函數(shù)平滑圖像，計(jì)算圖像的曲率和梯度.

(5)基于興趣區(qū)域設(shè)定初始輪廓和迭代次數(shù)，通過(guò)求解偏微分方程驅(qū)動(dòng)水平集的演化，得到最終的分割圖像.

偏微分方程PDE可以利用標(biāo)準(zhǔn)梯度下降方法最小化能量函數(shù).能量函數(shù)公式為：

(3)

算法流程圖如圖1所示.

圖1 本文算法流程圖Fig.1 Flow chart of the algorithm in this paper

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

4.1 實(shí)驗(yàn)過(guò)程

為了驗(yàn)證本文算法，實(shí)驗(yàn)環(huán)境為：處理器Intel(R)Xeon(R)E5-2620，CPU2.0GHz，內(nèi)存32GB，64位Windows Server2008 R2操作系統(tǒng).實(shí)驗(yàn)從每幅圖像中隨機(jī)取樣100個(gè)像素點(diǎn)，Li方法中的參數(shù)為：δ=1.0, Δt=0.1,v=0.001*255*255,λ1=λ2=1.0,u=1.0.實(shí)驗(yàn)根據(jù)不同的數(shù)據(jù)集，采用了不同的迭代次數(shù)，實(shí)驗(yàn)1和實(shí)驗(yàn)3選取的圖像與實(shí)驗(yàn)2的圖像相比，噪聲強(qiáng)度相對(duì)較弱，迭代200次后，圖像邊界分割清晰，成功收斂到目標(biāo)區(qū)域.實(shí)驗(yàn)2的頭頂灰度圖，噪聲強(qiáng)，為避免邊緣泄露問(wèn)題，迭代選取400次.

實(shí)驗(yàn)1如圖2所示，圖像為103*101的血管圖，實(shí)驗(yàn)選取13 493個(gè)像素中的100個(gè)像素進(jìn)行近似估計(jì)，初始輪廓為感興趣的不規(guī)則區(qū)域，迭代次數(shù)為200次.圖2(a)為原始圖像，2(b)中紅色虛線為初始輪廓，2(c)中紅色實(shí)線為分割參考輪廓，2(d)為兩種算法迭代100次后的分割圖像，2(e)為兩種分割算法最終的分割結(jié)果.(d)、(e) 中黑色線為本文算法，白色線為參考輪廓，藍(lán)色線為基于最小化區(qū)域擴(kuò)展擬合能量的圖像分割算法.

實(shí)驗(yàn)2如圖3所示，圖像為128*128的頭頂灰度圖，實(shí)驗(yàn)選取16384個(gè)像素中的100個(gè)像素近似估計(jì)，初始輪廓為感興趣的不規(guī)則區(qū)域，迭代次數(shù)為400次.圖3(a)為原始圖像，3(b)中紅色虛線為初始輪廓，3(c)中紅色實(shí)線為分割參考輪廓，3(d)為兩種算法迭代200次后的分割圖像，3(e)為兩種分割算法最終的分割結(jié)果.3(d)、3(e)中黑色線為本文算法，白色線為參考輪廓，藍(lán)色線為基于最小化區(qū)域擴(kuò)展擬合能量的圖像分割算法.

實(shí)驗(yàn)3如圖4所示，圖像為128*128的酵母熒光顯微圖，實(shí)驗(yàn)選取16 384個(gè)像素中的100個(gè)像素近似估計(jì)，初始輪廓為感興趣的不規(guī)則區(qū)域，迭代次數(shù)為200次.圖4(a)為原始圖像，4(b)中紅色虛線為初始輪廓，4(c)中紅色實(shí)線為分割參考輪廓，4(d)為兩種算法迭代100次后的分割圖像，4(e)為兩種分割算法最終的分割結(jié)果.4(d)、4(e)中黑色線為本文算法，白色線為參考輪廓，藍(lán)色線為基于最小化區(qū)域擴(kuò)展擬合能量的圖像分割算法.(注：圖2-4原圖是彩色，讀者如需了解詳情，請(qǐng)到本刊網(wǎng)站查看電子文稿)

(a)原始圖像 (b)初始輪廓 (c)參考輪廓 (d)100次迭代 (e)結(jié)果對(duì)比圖2 血管圖分割Fig.2 The bleod vessel segmentation

(a)原始圖像 (b)初始輪廓 (c)參考輪廓 (d)200次迭代 (e)結(jié)果對(duì)比圖3 頭頂灰度圖分割Fig.3 The perietal region segmentation

(a)原始圖像 (b)初始輪廓 (c)參考輪廓 (d)100次迭代 (e)結(jié)果對(duì)比圖4 酵母熒光顯微圖分割Fig.4 The yeast micrograph segmentation

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

從圖2(e)、圖3(e)、圖4(e)分割結(jié)果來(lái)看，本文算法能夠快速精確地收斂到目標(biāo)區(qū)域，克服了基于最小化區(qū)域擬合擴(kuò)展能量的圖像分割算法對(duì)于像素灰度相近的圖像分割效果不理想，易產(chǎn)生邊界泄露的問(wèn)題.

為了更好的分析實(shí)驗(yàn)結(jié)果，本文從算法分割時(shí)間和結(jié)果圖像與參考圖像的相似度Dice兩方面進(jìn)行評(píng)價(jià)，見(jiàn)表1、表2.算法分割時(shí)間以s為單位.

表1 分割時(shí)間比較
Tab.1 Comparision of the segmentation time s

圖像基于最小化區(qū)域擬合能量的圖像分割算法本文算法血管圖9．4214．787頭頂灰度圖22．3686．615酵母熒光顯微圖3．5052．758

表2 Dice相似度比較
Tab.2 Dice similarity comparison

圖像基于最小化區(qū)域擬合能量的圖像分割算法本文算法血管圖0．3260．946頭頂灰度圖0．8360．845酵母熒光顯微圖0．9000．910

5 結(jié)束語(yǔ)

本文針對(duì)現(xiàn)有基于最小化擴(kuò)展擬合能量的圖像分割方法，對(duì)于圖像像素灰度值相近時(shí)，分割效果不明顯，在弱邊緣處易產(chǎn)生邊界泄露的現(xiàn)象，提出了基于Nystrom方法的水平集醫(yī)學(xué)圖像分割算法.通過(guò)理論分析和相關(guān)實(shí)驗(yàn)，證明該方法能夠快速精確地收斂到目標(biāo)區(qū)域，在相同的迭代次數(shù)中，本文方法能更快分割圖像，并能保證較高的相似度.

[1]錢(qián)蕓,張英杰.水平集的圖像分割方法綜述[J]. 中國(guó)圖象圖形學(xué)報(bào),2008(1):7-13.

[2]CASELLES V, KIMMEL R, SAPIRO G. Geodesic Active Contours[C]// International Conference on Computer Vision. IEEE Computer Society, 1997:694.

[3]吳杰,朱家明,陳靜. 基于模糊聚類水平集的醫(yī)學(xué)圖像分割方法[J]. 計(jì)算機(jī)科學(xué), 2015(S2):155-159.

[4]劉雅婧,宋余慶,廖定安,等. 基于核模糊聚類的變分水平集醫(yī)學(xué)圖像分割方法[J]. 計(jì)算機(jī)應(yīng)用研究,2013, 30(11):3 510-3 513.

[5]郭元丁. 基于圖論的心臟CT圖像分割的研究[D]. 杭州：浙江大學(xué),2016.

[6]CARBALLIDO-GAMIO J, BELONGIE S J, MAJUMDAR S. Normalized cuts for spinal MRI segmentation[C]// CARS 2002 Computer Assisted Radiology and Surgery. Springer Berlin Heidelberg, 2002:1 054-1 054.

[7]劉鎖蘭,王江濤,王建國(guó),等. 一種新的基于圖論聚類的分割算法[J]. 計(jì)算機(jī)科學(xué),2008,35(9):245-247.

[8]閆成新,桑農(nóng),張?zhí)煨? 基于圖論的圖像分割研究進(jìn)展[J]. 計(jì)算機(jī)工程與應(yīng)用, 2006, 42(5):11-14.

[9]張洪剛,陳光,郭軍. 圖像處理與識(shí)別[M]. 北京：人民郵電出版社,2006.

[10]印世樂(lè),曾志勇. 一種改進(jìn)的Nystrom譜聚類圖像分割算法[J]. 計(jì)算機(jī)與現(xiàn)代化,2014(4):20-23.

[11]FOWLKES C, BELONGIE S, FAN C, et al. Spectral grouping using the Nystrom method[J]. Pattern Analysis & Machine Intelligence IEEE Transactions on, 2004, 26(2):214-225.

[12]陳應(yīng)良,王士同,錢(qián)蓉. 基于Nystrom方法的圖像譜分割算法的聚類改進(jìn)[J]. 計(jì)算機(jī)工程與設(shè)計(jì),2008,29(13):3 399-3 401.

[13]LI C, KAO C Y, GORE J C, et al. Minimization of region-scalable fitting energy for image segmentation[J]. IEEE Transactionson Image Processing, 2008, 17(10):1 940-1 949.