高亞明 劉兆邦 陳 斌 李 銘 黃來劍

1(揚(yáng)州大學(xué)信息工程學(xué)院 江蘇 揚(yáng)州 225127)2(中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所 江蘇 蘇州 215163)3(溫州市人民醫(yī)院 浙江 溫州 325699)

0 引 言

近年來，我國人民的生活質(zhì)量逐步提高，尿路結(jié)石的發(fā)病率也隨之遞增，南方發(fā)病率高于北方，約為22%～45%，個(gè)別省市甚至高達(dá)50%以上，且25%的患者需要住院治療。顯而易見，尿路結(jié)石已成為重要的公共衛(wèi)生問題[1]。國內(nèi)大多數(shù)研究表明，尿路結(jié)石的形成與很多因素有關(guān)，如地理環(huán)境、氣候、飲食習(xí)慣等。本文主要研究體內(nèi)結(jié)石中的單純性結(jié)石，其成分以草酸鈣最為多見，占結(jié)石的80%以上，無水尿酸次之[2]。目前利用傳統(tǒng)的結(jié)石成分分析方法對(duì)兩者進(jìn)行區(qū)分需要借助儀器，步驟復(fù)雜且成本較高，無法實(shí)現(xiàn)術(shù)前體內(nèi)無創(chuàng)檢測(cè)。因此，本文通過CT圖像，結(jié)合機(jī)器學(xué)習(xí)算法對(duì)兩者的結(jié)構(gòu)進(jìn)行分析鑒別[3]。

就目前而言，基于機(jī)器學(xué)習(xí)的輔助診斷算法已經(jīng)從純粹的理論發(fā)展到了臨床試驗(yàn)，在諸多方面如乳腺、肺結(jié)節(jié)等診斷上證實(shí)了其可行性與適用性[4]，但其在結(jié)石成分鑒別上的研究較少，比如Perrot等[5]利用放射和機(jī)器學(xué)習(xí)能準(zhǔn)確鑒別腎結(jié)石和靜脈增生，卻無法有效地在結(jié)石間進(jìn)行鑒別。因此，針對(duì)體內(nèi)單純性結(jié)石特征，本文將輔助診斷算法用于其成分分析中，輔助醫(yī)生進(jìn)行診斷。

1 圖像預(yù)處理

預(yù)處理主要運(yùn)用圖像增強(qiáng)技術(shù)，傳統(tǒng)的圖像增強(qiáng)技術(shù)大多是基于空間域?qū)D像進(jìn)行處理，本文利用灰度調(diào)整和圖像插值方法。

1.1 灰度調(diào)整

灰度調(diào)整利用灰度直方圖得到圖像中的像素亮度分布情況，再通過均衡化、規(guī)范化處理，調(diào)整圖像的清晰度。

本文首先將CT圖像的普通像素值更改為醫(yī)學(xué)專用的CT值，再根據(jù)結(jié)石組織CT值范圍截取[100,1 500]HU以內(nèi)的內(nèi)容進(jìn)行灰度調(diào)整[6]，其中HU(Hounsfiled Unit)值表示組織對(duì)X射線的吸收程度。調(diào)整效果如圖1所示。

圖1 灰度調(diào)整

1.2 圖像間插值

圖像間插值在每幀圖像中增加插值圖像以縮小距離，并將插值圖像和原始圖像共同組成三維數(shù)據(jù)以提高三維圖像的質(zhì)量[7]。

本文利用立方插值算法在圖像序列間進(jìn)行超分辨率重建，將CT圖像的像素間距統(tǒng)一為0.7 mm×0.7 mm×0.7 mm，如圖2所示。

圖2 圖像間插值

2 圖像分割

當(dāng)前醫(yī)學(xué)圖像分割正從手動(dòng)或半自動(dòng)向全自動(dòng)發(fā)展[8]。本文將臨床專家人工分割的結(jié)果作為金標(biāo)準(zhǔn)，利用3D Slicer軟件對(duì)病灶組織手動(dòng)勾勒感興趣區(qū)域(ROI)[9]并對(duì)其裁剪得到更精準(zhǔn)的ROI窗口[10]，分割過程如圖3所示。

圖3 圖像分割

3 特征提取

由于體外的草酸鈣結(jié)石質(zhì)硬，呈環(huán)形或桑葚形，表面粗糙有刺，而無水尿酸結(jié)石質(zhì)硬且表面光滑，呈圓形或卵圓形[11]，故本文針對(duì)體內(nèi)結(jié)石CT圖像，分別提取在二維和三維上的灰度、形狀和紋理特征[12]，并對(duì)其進(jìn)行對(duì)比分析。

3.1 灰度特征

醫(yī)學(xué)圖像中的灰度直方圖考慮形成圖像像素的強(qiáng)度，定義如下：由f(x,y)表示一幅灰度圖像，(x,y)處的值表示該位置像素，一幅圖有M個(gè)像素，分布在0～L-1灰度級(jí)，灰度直方圖則統(tǒng)計(jì)每級(jí)灰度像素?cái)?shù)目得到的統(tǒng)計(jì)圖。

(1)

式中：i表示灰度級(jí);L表示灰度級(jí)種類數(shù)(L≤256);mi表示灰度級(jí)為i的像素個(gè)數(shù)；M表示圖像總像素個(gè)數(shù)。

本文利用灰度直方圖，得到的灰度特征如下：

1)最大灰度值：分布在CT圖像中的最大灰度級(jí)。

2)均值：圖像灰度平均值。

(2)

3)方差：圖像灰度在數(shù)值上的分布情況。

(3)

4)熵：直方圖分布的均勻性。

(4)

5)傾斜度：直方圖分布的不對(duì)稱程度。

(5)

6)峰度：圖像的灰度分布在接近均值時(shí)的大致狀態(tài)。

(6)

3.2 形狀特征

形狀特征為結(jié)石分析提供重要的結(jié)構(gòu)信息，分為輪廓特征和區(qū)域特征。輪廓特征主要針對(duì)結(jié)石的外邊界，本文得到基于輪廓的特征如下：

1)緊密度：衡量一個(gè)形狀的緊致程度。

(7)

式中：A為形狀面積；i1和i2為二階矩：

式中：image(i,j)為CT圖像；(i,j)處的值為CT圖像中該位置的像素值；I為CT圖像中所有像素點(diǎn)的橫坐標(biāo)集合；J為CT圖像中所有像素點(diǎn)的縱坐標(biāo)集合。

2)最長(zhǎng)徑：邊界上相距最遠(yuǎn)的兩個(gè)點(diǎn)之間的距離。

3)離心率：用焦點(diǎn)間的距離除以長(zhǎng)軸的長(zhǎng)度。

區(qū)域特征則針對(duì)整個(gè)結(jié)石區(qū)域，本文得到基于區(qū)域的特征包括二維區(qū)域中的面積和三維區(qū)域中的體積。

3.3 紋理特征

紋理特征反映結(jié)石表面的粗糙度、光滑性[13]，提取方法分為統(tǒng)計(jì)分析、模型分析、結(jié)構(gòu)分析和頻譜分析。

本文采用基于統(tǒng)計(jì)的方法，主要包括：

1)灰度共生矩陣(GLCM)：統(tǒng)計(jì)圖像中不同灰度值的像素對(duì)出現(xiàn)的情況。

2)灰度游程矩陣(GLRLM)：統(tǒng)計(jì)圖像中具有相同灰度的直線區(qū)域(灰度、方向、長(zhǎng)度)出現(xiàn)的情況。

3)灰度區(qū)域尺寸矩陣(GLSZM)：統(tǒng)計(jì)圖像中具有相同灰度的區(qū)域(灰度、大小)出現(xiàn)的情況。

4)鄰域灰度差分矩陣(NGTDM)：統(tǒng)計(jì)圖像中相鄰區(qū)域的灰度差異情況。

5)灰度差分統(tǒng)計(jì)(GRAY_DIFF)：反映目標(biāo)場(chǎng)景在連續(xù)時(shí)間點(diǎn)圖像相減所構(gòu)成的圖像特征。

其中灰度共生矩陣(GLCM)最為經(jīng)典，應(yīng)用最廣，由它提取出如下特征：

1)角二階矩(energy)，反映圖像區(qū)域的均勻性或平滑性，表示為：

(8)

2)對(duì)比度(contrast)，表示圖像的紋理清晰度，表示為：

(9)

3)相關(guān)系數(shù)(correlation)，反映矩陣行與列的線性相關(guān)程度，值越大圖像區(qū)域灰度分布越均勻，表示為：

(10)

式中：μx、μy、σx、σy分別定義為：

4)熵(entrop)，度量圖像內(nèi)容的隨機(jī)性，表示為：

(11)

5)逆差分矩(uniformity)，定義為：

(12)

6)同質(zhì)化(homogeneity)，定義為：

(13)

7)和平均(sum average)，定義為：

(14)

式中：k=i+j。

所提特征如表1所示。

表1 所提特征

特征提取偽代碼如下：

功能：特征提取。

輸入：路徑字符串1，路徑字符串2。

輸出：2D灰度特征，2D形狀特征，2D紋理特征，3D灰度特征，3D形狀特征，3D紋理特征。

從磁盤讀取dcm醫(yī)學(xué)圖像數(shù)據(jù)文件

從磁盤讀取nrrd醫(yī)學(xué)圖像標(biāo)簽文件

forx=3 topatient_num

獲取dcm文件與nrrd文件的數(shù)量img_num

fori=1 toimg_num

按照InstanceID對(duì)dcm文件進(jìn)行排序，保證與nrrd文件順序一致

end

fori=1 toimg_num

查看當(dāng)前幀是否存在標(biāo)簽

ifm1>0 then

對(duì)當(dāng)前幀圖像預(yù)處理

fory=1 tonum

對(duì)當(dāng)前幀圖像的每個(gè)連通區(qū)域進(jìn)行分割，得到ROI

對(duì)ROI進(jìn)行裁剪得到最大ROI，即CT

將所有圖像Im和標(biāo)簽mm分別保存為三維矩陣O_CT3V、MSK3V

end

end

end

對(duì)CT提取二維灰度、形狀、紋理特征

將所提特征歸類為二維全局特征和內(nèi)部特征

對(duì)O_CT3V，MSK3V進(jìn)行裁剪

計(jì)算維度dim

ifdim>1

對(duì)三維圖像預(yù)處理

fory=1 tonum

對(duì)每個(gè)三維連通區(qū)域進(jìn)行分割，得到ROI

遍歷得到最大的ROI，并存儲(chǔ)為CT3V

end

forz=1 tosize(CT3V,3)

過濾，剔除沒有標(biāo)簽的幀圖像

end

對(duì)CT3V提取三維灰度、形狀、紋理特征

將所提特征歸類為三維全局特征和內(nèi)部特征

else

只有一幀圖像，直接將三維當(dāng)作二維處理，得到特征

end

end

4 特征選擇

特征選擇是對(duì)高維特征進(jìn)行篩選以降低特征維度的過程[14]，由于最小冗余最大相關(guān)(mRMR)算法綜合考慮了相關(guān)性和冗余性，故本文使用mRMR算法。

mRMR算法通過計(jì)算特征之間和類標(biāo)簽之間的互信息來選出冗余性最小和相關(guān)性最大的前N個(gè)特征，即從最小冗余和最大相關(guān)出發(fā)，給出一種基于互信息的評(píng)價(jià)準(zhǔn)則。其中：最大相關(guān)指的是特征與類標(biāo)簽信息之間的相關(guān)度最大；最小冗余指的是特征之間的冗余度最小[15]。

此外，將mRMR算法與相關(guān)系數(shù)(Spearma)法、Relief算法、SVM_RFE算法進(jìn)行對(duì)比分析，從而驗(yàn)證其在本文中的適用性。

5 分 類

本文采集的CT圖像數(shù)據(jù)是有限的，而支持向量機(jī)(SVM)分類器本身能較好地解決小樣本、高維度等問題，故本文選擇SVM分類器，針對(duì)有限樣本進(jìn)行訓(xùn)練。

在核函數(shù)的選擇上，基于高斯徑向基核函數(shù)的高維映射能力強(qiáng)、樣本適用性好、多用于特征數(shù)量小(15維)、樣本數(shù)量正常(119例)的情況，故本文使用基于高斯徑向基核函數(shù)的支持向量機(jī)(RBF_SVM)分類器。具體實(shí)現(xiàn)則采用第三方LIBSVM工具包，利用其內(nèi)置參數(shù)進(jìn)行調(diào)參。

同時(shí)，將其與基于Linear核函數(shù)的SVM(Linear_SVM)分類器，隨機(jī)森林(RF)分類器和Adaboost分類器進(jìn)行對(duì)比分析，從而驗(yàn)證RBF_SVM分類器在本文中的適用性。

6 實(shí) 驗(yàn)

6.1 實(shí)驗(yàn)數(shù)據(jù)和環(huán)境

1)實(shí)驗(yàn)所用數(shù)據(jù)：

從溫州市人民醫(yī)院采集結(jié)石CT圖像DICOM數(shù)據(jù)集，其中：草酸鈣患者59例；無水尿酸患者60例；層內(nèi)像素分辨率為512×512，層距為5 mm。

2)實(shí)驗(yàn)所用機(jī)器環(huán)境：

操作系統(tǒng)：Windows 10；CPU：Intel(R)Core(TM)i5-7500 CPU @ 3.40 GHz(3 401 MHz)；內(nèi)存：16 GB；編程軟件：MATLAB R2018a。

6.2 評(píng)估方法

1)分層10次10折交叉驗(yàn)證法。分層10折交叉驗(yàn)證首先將數(shù)據(jù)集分成10份，每份中類別之間的比例與整個(gè)數(shù)據(jù)集中的比例相同，輪番將其中9份作為訓(xùn)練集，1份作為測(cè)試集進(jìn)行實(shí)驗(yàn)，每輪得到一個(gè)精度，最后對(duì)10輪結(jié)果取平均作為最后的模型精度。模型精度估計(jì)需要進(jìn)行多次10折交叉驗(yàn)證，本文進(jìn)行10次10折交叉驗(yàn)證，對(duì)得到的10個(gè)精度(ACC)求均值。

2)調(diào)參。參數(shù)設(shè)定直接影響算法性能，本文利用交叉驗(yàn)證法，結(jié)合ROC曲線下的面積(AUC)對(duì)各分類器參數(shù)設(shè)置步長(zhǎng)，進(jìn)行調(diào)參，以選取最優(yōu)參數(shù)，優(yōu)化分類模型。

3)性能度量?；煜仃囀且粋€(gè)誤差矩陣，由它可得許多度量指標(biāo)，如表2所示。

表2 混淆矩陣

(15)

(16)

(17)

(18)

(19)

(20)

ROC曲線的橫軸為假正例率：

(21)

ROC曲線的縱軸為真正例率：

(22)

根據(jù)分類器的預(yù)測(cè)概率迭代更新閾值，得到多個(gè)混淆矩陣對(duì)應(yīng)的坐標(biāo)點(diǎn)，依據(jù)坐標(biāo)點(diǎn)畫出ROC曲線進(jìn)而計(jì)算AUC值。

6.3 實(shí)驗(yàn)步驟與結(jié)果

如圖4所示，首先對(duì)二維圖像進(jìn)行灰度調(diào)整，對(duì)三維圖像序列進(jìn)行插值，完成超分辨率重建，統(tǒng)一分辨率；其次，利用已標(biāo)好病灶的NRRD標(biāo)簽文件進(jìn)行圖像分割，并裁剪得到更精準(zhǔn)的ROI窗口；然后，再對(duì)ROI分別提取二維和三維上的灰度、形狀、紋理特征，并通過特征選擇算法進(jìn)行篩選；最后將篩選得到的特征放入分類器進(jìn)行模型訓(xùn)練，通過ACC、AUC等指標(biāo)對(duì)比分析分類效果。

圖4 實(shí)驗(yàn)步驟

1)特征提取算法對(duì)比。如表3和表4所示，各領(lǐng)域信息互補(bǔ)，綜合訓(xùn)練模型的鑒別性能，在統(tǒng)一使用mRMR算法與RBF_SVM分類器的前提下，紋理特征比灰度和形狀特征鑒別能力更強(qiáng)，2D特征略優(yōu)于3D特征，而提取所有特征時(shí)，ACC值和AUC值最高，即最能反映病灶組織空間分布的異質(zhì)性。

表3 特征提取平均ACC比較 %

表4 特征提取平均AUC比較 %

2)特征選擇算法對(duì)比。對(duì)各特征選擇算法進(jìn)行對(duì)比分析，由表5可知，在統(tǒng)一提取所有特征與使用RBF_SVM分類器的情況下，mRMR綜合考慮了特征之間的冗余性和特征類標(biāo)之間的相關(guān)性，故其選擇效果最好，驗(yàn)證了mRMR在本文中的適用性。

表5 特征選擇對(duì)比 %

3)分類器對(duì)比。由表6可知，在統(tǒng)一提取所有特征與使用mRMR的情況下，RBF_SVM因其在小樣本上的適用性較強(qiáng)而分類效果最好，驗(yàn)證RBF_SVM在本文中的適用性。

表6 分類器對(duì)比 %

4)最終輔助診斷算法評(píng)估。對(duì)結(jié)石CT圖像預(yù)處理及分割后提取所有特征，共130維，再利用mRMR算法篩選，篩選特征數(shù)對(duì)比如圖5所示。

圖5 篩選特征數(shù)對(duì)比

由圖5可知，選擇前15個(gè)特征分類效果最好，所選特征如表7所示。

表7 所選特征

將所選特征放入RBF_SVM分類器，得到本文所述的完整輔助診斷算法模型，綜合各項(xiàng)指標(biāo)對(duì)結(jié)果進(jìn)行評(píng)估，評(píng)估結(jié)果如表8所示。

表8 輔助診斷算法的各項(xiàng)指標(biāo) %

如圖6和圖7所示，10次結(jié)果的ACC值約為0.81，說明分類正確的樣本數(shù)較多，而ROC曲線對(duì)應(yīng)的AUC值在0.89左右，說明分類器性能較優(yōu)，所以，基于該輔助診斷算法可以實(shí)現(xiàn)對(duì)草酸鈣和無水尿酸結(jié)石的準(zhǔn)確分類。

圖6 輔助診斷算法的ACC,AUC指標(biāo)

圖7 輔助診斷算法的ROC曲線

7 結(jié) 語

本文結(jié)合臨床醫(yī)學(xué)，利用CT圖像進(jìn)行定量分析，并對(duì)比機(jī)器學(xué)習(xí)特征提取、選擇及分類算法的實(shí)驗(yàn)結(jié)果，提取結(jié)石在二維和三維上的灰度、形狀、紋理特征。最終使用mRMR特征選擇算法及RBF_SVM分類器，實(shí)現(xiàn)患者體內(nèi)草酸鈣和無水尿酸結(jié)石的準(zhǔn)確分類，ACC和AUC分別達(dá)到81.76%和89.03%，比現(xiàn)有結(jié)石分析方法更為快速簡(jiǎn)便，真正實(shí)現(xiàn)體內(nèi)無創(chuàng)，且符合臨床統(tǒng)計(jì)學(xué)結(jié)果，為臨床醫(yī)生診斷提供更為有效的參考依據(jù)。下一步將優(yōu)化圖像處理和機(jī)器學(xué)習(xí)算法，得到更高的準(zhǔn)確率。