血管粘連型肺結(jié)節(jié)圖像的序列分割方法

張 偉，張小龍，趙涓涓+，強 彥，唐笑先

(1.太原理工大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，山西 晉中 030600；2.賓夕法尼亞州立大學(xué) 信息科學(xué)與技術(shù)學(xué)院，賓西法尼亞州 尤尼弗西蒂帕克 16802；3.山西省人民醫(yī)院 PET/CT中心，山西 太原 030024)

0 引 言

肺癌是當(dāng)今發(fā)病率和死亡率最高的惡性腫瘤[1]。早期肺癌的5年生存率超過70%[2]。因此，早發(fā)現(xiàn)早治療就顯得尤為重要。而肺部病灶中，血管粘連型結(jié)節(jié)的惡性度非常大，其早期CT影像的病變區(qū)域卻難以被準確的分割，導(dǎo)致后續(xù)的識別和診斷出現(xiàn)誤差。同時，為了提高肺結(jié)節(jié)的檢出率，在CT掃描中常常采用超微小間隔CT掃描，因而又出現(xiàn)影像數(shù)據(jù)爆炸式增長與人工診斷力量嚴重不足的問題。因此，當(dāng)前迫切需要一種既高效又能準確分割該類型結(jié)節(jié)的方法。

1 相關(guān)工作

肺結(jié)節(jié)圖像的準確分割是后續(xù)特征提取和良惡性診斷的基礎(chǔ)。因此，近年來對肺結(jié)節(jié)圖像分割算法的研究日益受到人們重視。

Sun等[3]提出一種流向特征熵和測地線距離的K均值聚類算法來分割肺結(jié)節(jié)，取得較好的結(jié)節(jié)分割結(jié)果，但是隨著結(jié)節(jié)半徑增加，算法執(zhí)行效率就越低；Wei Y等[4]先利用多尺度Hessian矩陣對疑似結(jié)節(jié)進行增強，再通過設(shè)計mean-shift聚類的核函數(shù)，最終實現(xiàn)結(jié)節(jié)疑似感興趣區(qū)域(region of interest，ROI)的分割；Sivakumar S等[5]提出一種基于無監(jiān)督分割模型的possibilistic-fuzzy聚類算法，提高了肺結(jié)節(jié)分割的速度和準確率；鄒瑜等[6]針對自組織特征映射(self-organizing feature map，SOM)網(wǎng)絡(luò)隨著神經(jīng)元增加分割效果變差的問題，提出一種將有限脈沖響應(yīng)加入SOM中的方法，并通過合并聚類實現(xiàn)對乳腺超聲檢查圖像的有效分割；趙涓涓等[7]利用PET圖像中SUV閾值結(jié)合分水嶺分割實現(xiàn)了孤立性肺結(jié)節(jié)圖像的分割，分割效果較好。

雖然這些算法對實性結(jié)節(jié)的分割都很有效，但是血管粘連型結(jié)節(jié)的灰度值和血管的灰度值非常接近，分割結(jié)果中不能有效將血管和結(jié)節(jié)分離開。

因此，本文針對血管粘連型結(jié)節(jié)的序列CT圖像，提出了一種序列分割方法。與文獻中提到的其它方法進行比較，驗證了本方法的有效性和通用性。

2 CT圖像預(yù)處理

2.1 序列肺實質(zhì)圖像分割

該部分使用此前項目組提出的序列肺實質(zhì)分割方法[8]對CT序列圖像進行肺實質(zhì)分割。所使用的方法分為4個階段：第一，利用肺實質(zhì)圖像的位置特征得到肺部CT序列ROI圖像；第二，利用改進的超像素序列圖像分割算法對序列ROI圖像進行超像素分割；第三，采用遺傳算法(genetic algorithm，GA)優(yōu)化的自生成神經(jīng)森林(self-generating neural forest，SGNF)算法對超像素樣本進行聚類；第四，根據(jù)聚類后超像素集的灰度特征和位置特征識別肺實質(zhì)。

2.2 感興趣區(qū)域提取

該部分采用Song J等[9]提出的改進的雪橇算法實現(xiàn)肺結(jié)節(jié)的準確定位，并以此結(jié)節(jié)為中心，提取30×30的感興趣區(qū)域ROI，排除較大部分背景噪聲干擾的同時顯著提高后續(xù)圖像分割的效率。

3 改進的超像素序列分割

超像素[10]是由一系列位置相鄰且顏色、亮度、紋理等特征相似的像素點組成的小區(qū)域。其把具有相似特征的像素集合起來進行操作，降低了圖像處理的復(fù)雜度。

本文充分考慮CT圖像中血管和結(jié)節(jié)類型的特征，提出自適應(yīng)相似度系數(shù)的超像素序列分割算法(adaptive similarity coefficient sequential linear iterative clustering，ACSLIC)，算法使用Lab顏色空間中的L，a，b這3種顏色信息和像素點的x，y坐標構(gòu)成5維特征向量[L，a，b，x，y]T來表示CT圖像的單個像素點，并在像素間距離計算中引入自適應(yīng)相似度系數(shù)，設(shè)計新的適用于肺部CT圖像像素間相似性衡量的距離計算方法，如式(1)～式(3)所示

(1)

(2)

ds=w*dLab+dxy

(3)

其中，ds表示系數(shù)疊加后CT圖像兩個像素點間的相似度；dLab表示像素與聚類中心間的顏色特征距離；dxy表示像素與聚類中心間的空間位置特征距離；k是超像素的聚類中心，i是其搜索范圍內(nèi)某一像素點；w是k、i像素點間的自適應(yīng)相似度系數(shù)，其由式(8)所確定。

在CT圖像中，假設(shè)λ1和λ2為Hessian矩陣的兩個特征值，則線性結(jié)構(gòu)滿足條件λ1≈0，λ2?0，圓形結(jié)構(gòu)滿足λ1≈λ2?0。并且依據(jù)肺結(jié)節(jié)的類圓形特征和尺度不一特征。本文將Hessian矩陣與高斯函數(shù)相結(jié)合，通過設(shè)計多尺度的結(jié)節(jié)相似度函數(shù)來對CT圖像中的所有像素進行結(jié)節(jié)相似度值的計算，從而構(gòu)造自適應(yīng)的相似度系數(shù)w。其很大程度上拉大了血管和結(jié)節(jié)區(qū)域的相似度，有效避免了血管和結(jié)節(jié)被分割到同一個超像素樣本中，極大提高后續(xù)肺結(jié)節(jié)分割的準確率。如式(4)～式(8)所示

(4)

(5)

(6)

(7)

(8)

式(4)為2D圖像的Hessian矩陣定義；式(5)和式(6)實現(xiàn)圓形結(jié)構(gòu)的多尺度融合，其中σmin=dmin/4，σmin=dmax/4，dmin為肺結(jié)節(jié)圖像的最小直徑，dmax為肺結(jié)節(jié)圖像的最大直徑，實驗中設(shè)置dmin=3和dmax=30，N取8，并通過不斷的迭代尺度因子σ，選取各個尺度上響應(yīng)最大的點，此時具有最好的效果；式(7)為構(gòu)造的結(jié)節(jié)相似度函數(shù)，w1,w2表示R1和R2的權(quán)重系數(shù)，實驗中設(shè)置w1=0.5，w2=15；式(8)中No(λ)i表示第i個像素點的結(jié)節(jié)相似度值。

提出的ACSLIC算法的優(yōu)勢同樣還體現(xiàn)在如下兩個方面：

(1)采用近似符合肺結(jié)節(jié)類圓形形狀特征的六邊形進行聚類。在像素聚類的過程中不會以為不同方向上距離不同而導(dǎo)致聚類效果有所差異，從而分割得到的超像素樣本可以更好保留圖像的邊緣信息；

(2)精確定位肺結(jié)節(jié)，直接提取ROI進行操作，加快序列圖像分割速度的同時消除了較大部分背景噪聲的干擾；算法具體步驟如下：

算法1：改進的超像素序列分割算法(ACSLIC)

輸入：肺部ROI圖像序列；

輸出：超像素分割后的序列圖像；

(1)初始化六邊形網(wǎng)格中的K個聚類中心Ck=[Lk,ak,bk,xk,yk]T，六邊形網(wǎng)格間距S；

(2)執(zhí)行式(4)～式(8)，構(gòu)造顏色特征距離和空間位置特征距離間的自適應(yīng)相似度系數(shù)；

(3)在2*S為邊長的區(qū)域中，以[xk,yk]為開始點采用新構(gòu)造的相似度度量函數(shù)來搜尋與之相似的所有像素點，得到K個超像素樣本；

(4)計算所有超像素樣本中L、a、b、x、y的均值來更新聚類中心；

(5)重復(fù)(3)～(4)步，直到殘差小于固定的閾值T;

(6)對分割結(jié)果計算所有超像素樣本的鄰接矩陣，用于后續(xù)特征提取；

(7)輸入下一張待分割的ROI圖像，重復(fù)執(zhí)行(2)～(6)步；

(8)重復(fù)(1)～(7)步，直到ROI圖像序列分割完畢。

4 超像素樣本特征提取

4.1 提取新特征

4.1.1 對比度增強直方圖特征

灰度特征是肺結(jié)節(jié)區(qū)域與其余肺部組織最明顯的差異。然而肺部CT圖像通常因為噪聲、偽影等導(dǎo)致成像模糊、對比度低等問題，有必要在提取灰度特征之前對ROI圖像進行增強來突出像素間的差別,但是在圖像增強的過程中，會導(dǎo)致邊緣信息有所丟失。

因此，本文提出一種線性組合canny算子邊緣檢測圖像與直方圖均衡化圖像的方法，在增強ROI圖像對比度的同時保留了較多的圖像細節(jié)信息。具體步驟如下所示。

(1)計算ROI圖像各灰度級出現(xiàn)的概率，計算如式(9)所示

(9)

其中，nk表示ROI圖像中灰度級為gk的像素點的數(shù)目，L為最大的灰度級數(shù)。

(2)通過如式(10)的變換函數(shù)，建立gk到Gk的映射關(guān)系

(10)

(3)將原ROI圖像，標記為I；變換后的灰度值Gk對應(yīng)的圖像，標記為J1；

(4)利用canny算子提取ROI圖像的邊緣檢測圖像，標記為J2；

(5)將圖像J1和J2進行線性組合，作為最終圖像，標記為J；

(6)對圖像J進行超像素分割，提取每個超像素樣本對應(yīng)的灰度直方圖特征。

本文提出的方法與傳統(tǒng)的直方圖均衡化的對比如圖1所示。

圖1 本文方法與傳統(tǒng)直方圖均衡化的對比

4.1.2 超像素樣本鄰域紋理特征

對于原始ROI圖像，為了更好地辨別超像素樣本所處肺部組織的具體類型，比如結(jié)節(jié)類型、血管類型及背景類型。本文從突出每個超像素樣本與其鄰域超像素樣本的差異角度出發(fā)，設(shè)計了一組新特征來量化每個超像素樣本與其鄰接超像素樣本的紋理差異。

但是，由于：①血管粘連型結(jié)節(jié)的灰度值和血管的灰度值非常接近，不能有效將血管和結(jié)節(jié)分離開；②血管和結(jié)節(jié)區(qū)域的特征盡管很相似，但是在結(jié)節(jié)內(nèi)部往往會包含其它構(gòu)成成分(包含5種情況，分別為鈣化，軟組織，液體，脂肪和氣體)，并且結(jié)節(jié)在結(jié)構(gòu)上是類圓形，血管更趨向于線性結(jié)構(gòu)。因此，提取有效的紋理特征對區(qū)分血管和結(jié)節(jié)具有很大用處；③但是單個超像素的紋理特征不能有效提取血管和結(jié)節(jié)的紋理特征，這是因為血管和結(jié)節(jié)的紋理變化通常來自比超像素更大面積的區(qū)域，類似于文獻[11]中闡述。

因此，本文提出一種超像素樣本鄰域紋理特征，具體步驟如下：

圖2 高斯金字塔圖像與差異圖像的計算

(2)由于Level間的尺度不同，因此，先通過雙線性插值調(diào)整Level2，使其與Level1尺度相同；

(3)計算Level1與Level2的差異圖像，標記為S1，如圖2(c)(S1)所示；

(4)類似步驟(2)和步驟(3)，計算Level1與Level3、Level2與Level3以及Level2與Level4間的差異圖像，分別標記為S2，S3，S4。

(5)調(diào)整S1-S4為原ROI圖像大小，并對其分別執(zhí)行改進后的超像素分割算法。對每個超像素樣本分別提取灰度均值與灰度方差特征。計算如式(11)和式(12)所示

(11)

(12)

其中，ui(j)為圖像Sj中第i個超像素樣本的灰度均值；vari(j)為圖像Sj中第i個超像素樣本的灰度方差；因此，實驗中標記單個圖像中單個超像素樣本對應(yīng)的特征為SFi=[ui,vari]。

(6)構(gòu)造當(dāng)前超像素樣本及其鄰近左上、左下、右下、右上4個超像素樣本的特征集合，標記為ASFi=[SFi,SFi1,SFi2,SFi3,SFi4]，如圖3所示。從而構(gòu)成最終4×2×5=40維的特征向量。

圖3 超像素樣本的鄰域關(guān)系

4.1.3 位置信息特征

除辛北斷層外，研究區(qū)內(nèi)的斷層級別稍低，延伸也較短，在沙三下亞段至沙四段形成了一些地塹和地壘。除X176西的南北向斷層外，這些斷層大都向下貫穿至沙四下-孔店組，雖然這些斷層活動強度均較弱，斷距10～50m，對沙四段沉積無明顯的控制作用，但由于其切穿沙四下亞段的鹽膏層，在斷層活動期，仍可對流體起到斷裂輸導(dǎo)作用。

本文根據(jù)提取的ROI圖像的中心坐標必處于肺結(jié)節(jié)區(qū)域，從而提出一種新的超像素樣本位置信息特征。表示如式(13)和圖4所示

(13)

其中，(xc,yc)為ROI圖像的中心坐標；(xi,yi)為第i個超像素樣本的中心坐標；H和W分別為ROI圖像的高度和寬度。

圖4 ROI圖像中心與超像素樣本中心的距離關(guān)系

4.2 特征數(shù)據(jù)歸一化

特征數(shù)據(jù)的歸一化是后續(xù)聚類的預(yù)處理步驟。歸一化的主要目的是：有效避免具有更大數(shù)值范圍特征的主導(dǎo),并在聚類的過程中顯著降低數(shù)值計算的困難,因此，本文采用如式(14)所示的方式對特征數(shù)據(jù)進行歸一化

(14)

其中，F(xiàn)是一個原始的特征值；F′是歸一化后的特征值；從而，特征向量中的每個值的取值范圍都為[0,1]。

5 距離約束稀疏子空間聚類

稀疏子空間聚類(sparse sub-space clustering，SSC)[12]是一種基于譜聚類的子空間聚類方法。給定一組數(shù)據(jù)，該聚類方法先將所有數(shù)據(jù)作為字典來線性表示每個數(shù)據(jù)，接著采用表示系數(shù)矩陣的稀疏約束使數(shù)據(jù)盡可能被同一子空間的數(shù)據(jù)線性表示，并構(gòu)造關(guān)于子空間系數(shù)的相似度矩陣，最后利用譜聚類等方法得到數(shù)據(jù)聚類的結(jié)果。

在本文中，提取的ROI圖像是以結(jié)節(jié)為中心的，并且肺結(jié)節(jié)圖像具有類圓形狀特性?；诖讼闰炛R，提出一種基于距離約束的稀疏子空間聚類方法。距離約束的引入有利于結(jié)節(jié)數(shù)據(jù)能夠更好的被同一子空間內(nèi)鄰近的結(jié)節(jié)特征數(shù)據(jù)線性表示，進一步有效排除血管對結(jié)節(jié)聚類的影響。采用的距離約束函數(shù)如式(15)所示

(15)

其中，Di和Dj為第i和第j個超像素樣本到ROI圖像中心的距離。引入距離約束的稀疏子空間聚類模型，如式(16)所示

(16)

(17)

其中，u為子空間表示系數(shù)；U=[u1,u2,…,uK]∈RK×K，為表示系數(shù)矩陣；距離越大，1/wji越小，否則1/wji越大；W為相似度矩陣；因此，距離約束有利于結(jié)節(jié)數(shù)據(jù)盡可能被同一子空間鄰近的最相似的數(shù)據(jù)線性表示。具體的算法描述如下：

算法2：距離約束稀疏子空間聚類

輸入：待聚類圖像I，分割后的超像素樣本總個數(shù)K；

輸出：聚類后的肺結(jié)節(jié)圖像；

(1)ACSLIC分割得到K個超像素樣本；

(2)對每個超像素樣本提取對比度增強直方圖特征、鄰域紋理特征以及位置信息特征。得到N個特征矩陣；

(3)利用式(15)計算超像素間的距離約束值；

(4)利用式(16)得到子空間表示系數(shù)矩陣U；

(5)利用式(17)構(gòu)造相似度矩陣W；

(6)根據(jù)W，采用歸一化割法(normalized cut，Ncut)算法得到超像素樣本聚類結(jié)果。

6 實驗結(jié)果與分析

本文算法的實驗平臺是Visual Studio 2013和MATLAB R2014b，PC處理器為Intel(R) Core(M) i7-4770，主頻3.40 GHz，內(nèi)存8 GB。實驗數(shù)據(jù)來源于LIDC-IDRI[13]和山西省人民醫(yī)院，挑選80組血管粘連型CT序列圖像共1468張，單張CT圖像大小為512×512。實驗中我們邀請到兩位影像科醫(yī)師做手動肺結(jié)節(jié)分割。同時為了驗證本文方法的有效性，針對所有的CT序列圖像，對本文方法、脈沖耦合神經(jīng)網(wǎng)絡(luò)算法(PCNN-pulse coupled neural network,PCNN)[14]、模糊c均值聚類算法(fuzzy c-means，F(xiàn)CM)[15]以及最新的流向特征熵和測地線距離的粘連血管型肺結(jié)節(jié)聚類分割方法(flowing entropy and geodesic distance，F(xiàn)EGD)[3]進行了大量的對比實驗。

6.1 定性分析

實驗中CT序列圖像數(shù)目較多,本文只對部分分割過程進行說明，如圖5所示。

同時，本文與FEGD、PCNN以及FCM算法進行了對比實驗。部分結(jié)果如圖6所示。

圖5 本文方法對血管粘連型結(jié)節(jié)圖像序列的分割結(jié)果注:(a)原始CT序列圖像;(b)分割得到的序列肺實質(zhì)圖像;(c)提取的ROIs圖像;(d)ACSLIC對ROIs圖像分割的結(jié)果;(e)距離約束稀疏子空間聚類得到的肺結(jié)節(jié)圖像掩膜序列;(f)本文方法的分割結(jié)果;(g)醫(yī)師手動分割的結(jié)果

圖6 4種方法對血管粘連型結(jié)節(jié)圖像分割結(jié)果的對比 注：(a)原始CT圖像；(b)對(a)局部放大的結(jié)果；(c)～(f)分別為PCNN、FCM、FEGD以及本文方法分割得到的肺結(jié)節(jié)圖像掩膜；(g)醫(yī)師手動分割的結(jié)果。

從結(jié)果圖可以直觀觀察到，PCNN和FCM對血管粘連型結(jié)節(jié)的分割結(jié)果不準確，不能有效分離血管和結(jié)節(jié)，導(dǎo)致分割得到的結(jié)節(jié)粘連有較多的血管組織；此外，F(xiàn)EGD的分割效果較好，結(jié)果優(yōu)于PCNN和FCM，但是分割結(jié)果存在少量的邊界泄露(如圖6，行1-4，行6所示)；并且FEGD在特定的CT圖像中會不可避免的引入部分血管(如圖6，行5所示)，這是因為FEGD不能有效處理血管橫斷面為圓形的血管；然而，本文方法可以有效的分離血管和結(jié)節(jié)，并且改進的超像素分割算法保留了更多的邊界信息，分割結(jié)果更平滑。相比這3種典型分割算法，本文方法都能取得最好的分割效果。并且與醫(yī)師手動分割的結(jié)果基本一致。

6.2 定量比較

本文通過采用變化信息、概率邊緣指數(shù)以及時間復(fù)雜度來對分割方法進行定量比較。

6.2.1 變化信息

變化信息(variation of information，VoI)是計算圖像像素點從一個聚類轉(zhuǎn)變?yōu)榱硗庖粋€聚類丟失或獲得信息量的大小[16]。假設(shè)參考分割圖像與實際分割結(jié)果分別用Sa和Sb表示。則可利用參考分割圖像Sa的熵[17]H(Sa)、實際分割結(jié)果Sb的熵H(Sb)以及Sa和Sb的聯(lián)合熵I(Sa,Sb)，得VoI計算如式(18)所示

VoI(Sa,Sb)=H(Sa)+H(Sb)-2I(Sa,Sb)

(18)

(19)

(20)

(21)

其中，P(m)表示一個像素被分割到第m個區(qū)域中的概率。VoI越小，表示實際分割結(jié)果與參考分割圖像越接近，分割效果越好。

6.2.2 概率邊緣指數(shù)

概率邊緣指數(shù)(probabilistic rand index，PRI) 檢驗實際分割結(jié)果與參考結(jié)果之間的屬性共生的一致性的參量[18]。假設(shè)(xi,yi)是原圖S中的一個像素對，如果在Sa中標記為(ai,aj)且相同，則在Sb中標記為(bi,bj)也應(yīng)該相同。PRI的計算如式(22)所示

(22)

其中，M為原始圖像S中的像素個數(shù)；I為判別函數(shù)，主要作用是判斷像素對是否具有相同標記。PRI 值越大，說明實際分割結(jié)果與醫(yī)師手動分割結(jié)果的屬性共生一致性越好，表明分割結(jié)果越好。

此外，為了驗證本文方法的通用性，實驗中計算了所有序列CT數(shù)據(jù)分割結(jié)果的VoI與PRI的均值。見表1。

表1 4種算法對所有血管粘連型結(jié)節(jié)圖像序列的VoI與RPI均值

實驗結(jié)果表明，對血管粘連型結(jié)節(jié)的CT序列圖像，PCNN和FCM兩種算法的分割效果基本一致；FEGD的分割效果遠優(yōu)于PCNN和FCM；但是，相比PCNN、FCM和FEGD這3種分割算法，本文方法都達到最低的VoI均值和最高的PRI均值，從而，從客觀的定量比較結(jié)果進一步驗證了本文方法的在分割這類型結(jié)節(jié)圖像上的優(yōu)越性。

6.2.3 時間復(fù)雜度

通過對所有的實驗數(shù)據(jù)進行分析，發(fā)現(xiàn)一個CT序列圖像包含結(jié)節(jié)的平均數(shù)量為8張。因此，表2顯示了4種算法序列分割的平均處理時間。雖然對于單個CT序列，平均處理時間比FCM算法多出0.48 s，但是本文方法的分割效果遠優(yōu)于FCM算法的分割結(jié)果，在醫(yī)學(xué)上這種用時間來換取分割的精度是值得的。

表2 4種算法分割序列CT圖像的平均處理時間

7 結(jié)束語

針對血管粘連型結(jié)節(jié)的序列CT圖像，本文提出一種基于超像素和稀疏子空間聚類的肺結(jié)節(jié)圖像序列分割方法，能夠?qū)崿F(xiàn)對血管粘連這種困難分割結(jié)節(jié)圖像的準確分割。該方法的主要貢獻是：在超像素分割算法中構(gòu)造了自適應(yīng)相似度系數(shù)；在超像素樣本特征提取中，設(shè)計提取了對比度增強直方圖特征、超像素樣本鄰域紋理特征以及位置信息特征等多個新特征；在稀疏子空間聚類算法中引入距離約束項；這3個方面都有效分離了血管和結(jié)節(jié)。從定性分析和定量比較結(jié)果可以看出，本文能更好地分離血管和結(jié)節(jié)，邊緣輪廓平滑，具有最好的序列分割結(jié)果，且與醫(yī)師手動分割結(jié)果基本一致。因此，本文的序列分割方法能更完整、準確地分割出血管粘連型結(jié)節(jié)圖像序列。