肝囊型包蟲病超聲圖影像區(qū)域分割算法研究

王正業(yè)，熱娜古麗·艾合麥提尼亞孜，王曉榮，米吾爾依提·海拉提，嚴(yán)傳波

1.新疆醫(yī)科大學(xué) a.公共衛(wèi)生學(xué)院；b.醫(yī)學(xué)工程技術(shù)學(xué)院，新疆 烏魯木齊 830011；2.新疆醫(yī)科大學(xué)第一附屬醫(yī)院 超聲診斷科，新疆 烏魯木齊 830011

引言

肝包蟲病作為一種人畜共患的寄生蟲病，在畜牧地區(qū)特別是醫(yī)療衛(wèi)生條件較差的偏遠(yuǎn)地區(qū)較為流行[1]。該疾病早期無明顯臨床癥狀，多數(shù)情況下患者在常規(guī)體檢時才被發(fā)現(xiàn)感染了包蟲病，而在醫(yī)療服務(wù)不發(fā)達(dá)的地區(qū)，大多數(shù)患者產(chǎn)生自覺癥狀時肝囊型包蟲病灶已發(fā)展至壓迫周圍組織器官的程度，這不僅增加了治療難度，同時還會使患者承受更多的痛苦，消耗額外的醫(yī)療資源。因此肝包蟲發(fā)病早期的精確診斷對于后續(xù)的臨床治療及預(yù)后有重要意義。超聲檢查作為一種簡單易行、無創(chuàng)、無輻射、精準(zhǔn)成像清晰的影像學(xué)方法，是肝囊型包蟲病診斷的首選檢查方法，其在偏遠(yuǎn)地區(qū)該疾病的早期篩查中發(fā)揮重要的作用[2]。由于醫(yī)師對于超聲圖像的診斷存在一定的主觀性，尤其是在偏遠(yuǎn)的畜牧地區(qū)，往往存在醫(yī)生經(jīng)驗不足的情況，影響肝包蟲發(fā)病早期的精確診斷，針對此類情況，有研究者提出使用計算機輔助診斷（Computer Aided Diagnosis，CAD）的方法輔助醫(yī)生診斷疾病，以此來彌補部分地區(qū)醫(yī)療資源不足的情況[3]。在肝包蟲病的CAD研究中，有研究者使用傳統(tǒng)機器學(xué)習(xí)方法和基于深度學(xué)習(xí)的方法在肝囊型包蟲病CT圖像病灶定位及類型區(qū)分任務(wù)中獲得了較好的結(jié)果[4-6]，但目前有關(guān)肝囊型包蟲病超聲圖像自動分割的研究較少。Wu等[7]使用手動分割的方法，提取肝包蟲病超聲影像病灶區(qū)域，并使用卷積神經(jīng)網(wǎng)絡(luò)（Convlutional Neural Networks，CNN）實現(xiàn)了肝囊型包蟲病超聲影像單囊型（CE-1）、多子囊型（CE-2）、內(nèi)囊塌陷型（CE-3）、壞死實變型（CE-4）、鈣化型（CE-5）5種分型的自動診斷，其中VGG19模型取得了最好的結(jié)果（精確度為90.6%）。鄒奕軒等[8]在甲狀腺結(jié)節(jié)良惡性自動診斷的研究中，手動去除了超聲圖像中除感興趣區(qū)域（Region of Interest，ROI）之外的信息，此外，許多有關(guān)超聲圖像的智能診斷研究均采用手動分割的方法來提升模型的性能[9-11]。人工標(biāo)注的方法雖然能較好地分割出目標(biāo)區(qū)域，但耗時較長，在實際應(yīng)用于邊緣畜牧區(qū)疾病篩查診療時較困難，且人工標(biāo)注截取ROI病灶區(qū)域時由于個體的主觀性差異可能會遺漏有利于模型分類的圖像信息，所以研究出一種自動分割肝囊型包蟲病超聲圖像（一般為扇形區(qū)域）的算法來提升計算機輔助診斷性能具有重要意義，此外在基于醫(yī)學(xué)影像的計算機輔助診斷系統(tǒng)中，超聲儀器在報告結(jié)果時生成的非影像信息會增大分類或病灶檢測模型的計算量，開展對肝囊型包蟲病超聲圖像扇形影像區(qū)域的自動分割算法研究，精準(zhǔn)的圖像分割對后續(xù)病灶分型有重要的意義，基于此，本文旨在使用Otsu閾值分割法[12]、馬爾可夫隨機場算法[13]及基于深度學(xué)習(xí)的Poly-YOLO算法[14]對囊型肝包蟲病超聲圖像中的扇形影像區(qū)域進(jìn)行分割，以期去除圖像中干擾信息，降低后續(xù)分類器訓(xùn)練時的計算量，進(jìn)而提升模型的整體性能。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)集采集

本研究選用的數(shù)據(jù)集來源于新疆醫(yī)科大學(xué)第一附屬醫(yī)院2008年1月至2018年12月采集的囊型肝包蟲病影像資料，本研究已通過新疆醫(yī)科大學(xué)第一附屬醫(yī)院倫理審查委員會批準(zhǔn)（K202207-12）。參考世界衛(wèi)生組織包蟲病非正式工作組關(guān)于包蟲病灶分型標(biāo)準(zhǔn)[15]，本次納入研究的肝囊型包蟲病超聲圖像分為5種類型：單囊型、多子囊型、內(nèi)囊塌陷型、壞死實變型、鈣化型。本研究收集了204例患者共計807張肝包蟲病超聲影像圖片，每張超聲圖像的類型都由該醫(yī)院擁有3年以上工作經(jīng)驗的超聲科影像醫(yī)師驗證確定。由于本研究納入的超聲圖片是由不同型號的儀器采集，在分析之前需要對圖像進(jìn)行預(yù)處理操作：① 尺寸歸一化：超聲儀器報告結(jié)果中除扇形成像區(qū)域外，往往會包含超聲儀器成像參數(shù)、品牌及醫(yī)院信息（圖1a），這些信息在基于影像的計算機輔助診斷過程中會影響判讀結(jié)果，增加計算量，影響模型的性能，因此將所有納入研究的圖像尺寸歸一化至840 px×720 px，便于后續(xù)分割算法的計算；② 圖像去噪：利用高斯濾波去除圖像中的小噪聲點，同時最大程度保留影像資料的邊緣信息，提升分割算法性能，并使用Python編程語言中的labelme圖像標(biāo)注軟件對圖片進(jìn)行標(biāo)注（圖1b），由于肝囊型包蟲病超聲影像的特殊性，特別是單囊型與多子囊型，其病灶區(qū)域不易與背景區(qū)分，病灶區(qū)域多為無回聲的液性暗區(qū)，僅有少量的回聲光點（圖1c）。

圖1 肝囊型包蟲超聲圖

1.2 方法

1.2.1 Otsu閾值分割法

Ostu算法是由日本學(xué)者Kumar等[16]提出的用于圖像分割領(lǐng)域的經(jīng)典算法，也被稱為最大類間方差算法，主要利用圖像的灰度分布確定一個合適的閾值K，將圖像分成前景和后景2個部分，即所需目標(biāo)（Objective，O）與背景（Background，B），并在O與B之間的類間方差取最大值，以此來獲得最合適的閾值，類間方差計算方式如公式（1）所示。

式中，μ代表圖像中所有像素灰度的平均值，μ0'與μb'代表目標(biāo)區(qū)域與背景區(qū)域的灰度均值，當(dāng)K的取值使得e2(K)最大時，即獲得了最合適的閾值K[17]。本研究采用單尺度Retinex（Single Scale Retinex，SSR）算法對超聲影像進(jìn)行增強處理，增強圖像前景與背景對比度，在此基礎(chǔ)上進(jìn)行閾值分割，以獲取較好的分割結(jié)果。SSR算法具體原理為將輸入的圖像F(x,y)分解為反射圖像R(x,y)和照射圖像L(x,y)，其中L(x,y)被認(rèn)為是影像圖像質(zhì)量的主要因素，R(x,y)被認(rèn)為是不受外界影響的原始圖像，其定義如公式（2）所示。

SSR算法即為在原圖F(x,y)的基礎(chǔ)上去除L(x,y)干擾因素，還原圖像的本質(zhì)信息，應(yīng)用于圖像處理領(lǐng)域時，即為對圖像的增強，表示方式如公式（3）所示。

式中，*為卷積算子，G(x,y)為標(biāo)準(zhǔn)環(huán)繞函數(shù)，通過處理F(x,y)中的像素與環(huán)繞加權(quán)均值之間的比值得到對于干擾信息L(x,y)的估計，最終得到R(x,y)的取值[18]。

1.2.2 馬爾可夫隨機場

馬爾可夫隨機場（Markov Random Field，MRF）模型是基于概率圖的數(shù)學(xué)模型，從聚類的角度實現(xiàn)對圖像的分割，其將具有相同性質(zhì)的像素點識別為同一類，即將分割問題轉(zhuǎn)化為標(biāo)簽分類問題，傳入圖像C={Ci}i∈S，位置集合為S，其中Ci代表圖像在i處的影像學(xué)特征，定義標(biāo)記集合 X={Xi}i∈S，對于位置 i處的標(biāo)記，Xi∈ {1,2…,k}，k為像素分類的類別數(shù)。根據(jù)貝葉斯（Bayesian）公式在標(biāo)記集合X中實現(xiàn)最優(yōu)影像分類的估計，對目標(biāo)函數(shù)進(jìn)行求解方式如公式（4）所示。

式中，P(X=x)為標(biāo)記集合的聯(lián)合概率分布，一般來說，目標(biāo)點的標(biāo)記值與領(lǐng)域像素點的標(biāo)記值是相關(guān)的，該聯(lián)合概率分布被認(rèn)為符合馬爾可夫模型[19]，P(C|X=x)代表傳入圖像C的條件分布概率，也被稱為似然函數(shù)，在一般情況下，似然函數(shù)可以被認(rèn)定為符合高斯分布。Hammersley-Clifford定理指出[20]，吉布斯隨機場（Gibbs）與馬爾可夫隨機場具有等價關(guān)系，此時確定Gibbs隨機場的能量函數(shù)就能夠確定馬爾可夫隨機場，求解方式如公式（5）～（8）所示。

式中，U（x）為能量函數(shù)，Z為歸一化常數(shù)，已知P(X=x)具有馬氏性，則進(jìn)一步推導(dǎo)如公式（6）所示。

式中，Ni表示位置i的領(lǐng)域集合。假定P(C|X=x)服從高斯分布，即各點之間相互獨立，則一步推導(dǎo)如公式（7）所示。

將公式（6）與（7）代入目標(biāo)函數(shù)，即可求得最優(yōu)解，求解方式如公式（8）所示。

1.2.3 Poly-YOLO算法

Poly-YOLO是在YOLOv3（You Look Only Once）[21]算法的基礎(chǔ)上改進(jìn)的一種基于深度學(xué)習(xí)的實例分割方法。與YOLOv3相比，Poly-YOLO消除了大量重寫標(biāo)簽與無效錨點分配的問題，同時使用邊界多邊形執(zhí)行實例分割，在訓(xùn)練及檢測與尺寸無關(guān)的多邊形后會生成具有不同數(shù)量頂點的多邊形，其結(jié)構(gòu)如圖2所示。在圖像特征提取層，采用了DarkNet-53的網(wǎng)絡(luò)[22]結(jié)構(gòu)，其采用了與Residual network結(jié)構(gòu)相似的殘差塊，不同之處在于特征提取的每層都使用了較少的卷積濾波器，并且采用了壓縮激發(fā)模塊（Squeeze and Excitation，SE）去擴(kuò)展它（圖 3），SE允許使用空間和通道的信息，使預(yù)測的準(zhǔn)確率進(jìn)一步提高。在網(wǎng)絡(luò)融合層采用輕量級的帶有Hypercolmn（HC）結(jié)構(gòu)的block替代，Hypercolmn作用為將原有的多尺度部分合成為單尺度輸出，此外使用Stairstep插值方法增加圖像的分辨率（圖4），與YOLOv3原有的直接上采樣相比，輸出更加平滑。

圖2 Poly-YOLO網(wǎng)絡(luò)結(jié)構(gòu)

圖3 改進(jìn)的殘差模塊

圖4 梯度上采樣結(jié)構(gòu)

Poly-YOLO訓(xùn)練環(huán)境：CPU為酷睿i7-7700K；GPU為GTX1080Ti；運行內(nèi)存16 GB；實驗用操作系統(tǒng)為ubuntu20.04；在tensor flow-GPU框架上進(jìn)行開發(fā)訓(xùn)練；Batchsize設(shè)置為8；最大迭代次數(shù)設(shè)置為200，當(dāng)loss值經(jīng)過10次迭代后不再下降，停止訓(xùn)練。

1.3 評價指標(biāo)

本研究采用Dice相似系數(shù)（Dice Similarity Coefficient，DSC）、重疊度（Intersection of Union，IOU）、真陽性率（True Positive Rate，TPR）、豪斯多夫距離（Hausdorff Distance，HD）評價上述3種算法。其中DSC為一種集合相似度量指標(biāo)，可用來計算2個樣本的相似度，取值范圍為0～1，該值越接近1代表兩樣本相似度越高，反之則相似度越低，計算方式如公式（9）所示；IOU分?jǐn)?shù)作為一種分割任務(wù)中廣泛應(yīng)用的標(biāo)準(zhǔn)性能度量指標(biāo)，用于評價目標(biāo)區(qū)域和真實區(qū)域之間的相似性，計算方式如公式（10）所示；TPR為預(yù)測正確的樣本占總陽性樣本的比例，又稱為靈敏度，計算方式如公式（11）所示；HD是衡量1個集合到另1個集合中最近點的最大距離，用于檢測2個集合的不匹配程度，該值越小代表2個集合越相似，由于從集合A到集合B和集合B到集合A的HD不同，所以定義二者直接距離的計算方式如公式（12）所示。

式中，P1代表分割預(yù)測的結(jié)果；T1代表真實結(jié)果。

式中，TP、FP、FN分別代表真陽性數(shù)、假陽性數(shù)和假陰性數(shù)。

式中，TP、FN分別代表真陽性數(shù)和假陰性數(shù)。

1.4 統(tǒng)計學(xué)分析

本研究采用Python 3.6.0語言構(gòu)建相關(guān)的分割算法及模型，使用Python中Scipy1.5.0及Numpy 1.18.5庫對DSC、IOU、TPR、HD等評價指標(biāo)進(jìn)行計算，使用箱式圖描述分割算法評價指標(biāo)的結(jié)果，并利用Matplotlib 3.2.2庫及GraphPad Prism 8.4.2進(jìn)行相關(guān)圖表的繪制。

2 結(jié)果

2.1 Otsu閾值分割法

肝囊型包蟲超聲影像原圖如圖5a所示，單尺度SSR圖像增強的算法處理超聲圖片如圖5b所示，在原圖與經(jīng)SSR處理的圖片上進(jìn)行Ostu閾值分割結(jié)果如圖5c～d所示。對于圖中微小連接處使用腐蝕算法進(jìn)行處理，如圖6a～b所示，采用最大連通域算法在閾值分割后的圖像上尋找面積最大閉合區(qū)域，再與原圖進(jìn)行邏輯“與”運算，得到最終分割結(jié)果（圖6c～e）。結(jié)果顯示，使用SSR算法對肝囊型包蟲病圖像處理發(fā)現(xiàn)目標(biāo)與背景的對比度明顯增強，使用Ostu分割后可以觀察到更明顯的差異，經(jīng)過增強的圖像能更完整地將影像區(qū)域進(jìn)行較好的分割，原圖對比度較低，Ostu分割方法無法較好地分割扇形區(qū)域；同時由于儀器報告格式的不同，部分超聲圖像存在影像區(qū)域與非影像區(qū)域連接現(xiàn)象。本文采用腐蝕算法對二值圖像進(jìn)行處理，將微小連接區(qū)域腐蝕，獲得完整且獨立的扇形。

圖5 SSR算法的應(yīng)用情況

圖6 Ostu閾值分割過程

2.2 馬爾可夫隨機場

本研究中，每張圖像分類類別設(shè)置為2種（目標(biāo)與背景），迭代次數(shù)設(shè)置為200次，經(jīng)過迭代，部分圖像出現(xiàn)目標(biāo)劃分為背景、背景劃分為目標(biāo)的情況（圖7a），對此類圖進(jìn)行像素反轉(zhuǎn)操作（圖7b），得到二值化圖像后，同樣利用最大連通域算法及邏輯與運算得到分割結(jié)果（圖7c～d）。本研究中，由于不同圖像之間像素灰度值分布有不同，導(dǎo)致2種類別的高斯密度函數(shù)存在差異，使得部分圖像的目標(biāo)與背景錯分，對此類圖片進(jìn)行灰度值反轉(zhuǎn)操作。雖然在該方法的研究中存在類別劃分誤差的問題，但由于馬爾可夫模型可以較好地將圖片中顏色、紋理相似部分歸為一類，適用于肝包蟲病超聲影像扇形區(qū)域分割任務(wù)，能夠較好地擬合超聲圖像扇形影像區(qū)域。得益于對目標(biāo)與背景的精確劃分，馬爾可夫模型在DSC和IOU性能評價指標(biāo)中具有良好的表現(xiàn)。

2.3 Poly-YOLO算法

Poly-YOLO經(jīng)過72輪epoch訓(xùn)練后終止訓(xùn)練，Loss值在62輪后不再下降，使用訓(xùn)練后的模型對肝囊型包蟲病超聲圖像原圖上進(jìn)行扇形區(qū)域預(yù)測，結(jié)果如圖8a所示，生成掩模圖（圖8b），通過與原圖與運算，得到扇形區(qū)域（圖8c）。與傳統(tǒng)方法不同，基于深度學(xué)習(xí)的Poly-YOLO算法需要大量的數(shù)據(jù)輸入到模型中進(jìn)行訓(xùn)練，使用極坐標(biāo)系來表述圖像中目標(biāo)區(qū)域的位置，學(xué)習(xí)極坐標(biāo)單元原點與目標(biāo)位置之間的距離，同時可以通過極坐標(biāo)單元內(nèi)的角度來調(diào)整位置。結(jié)果顯示分割區(qū)域并不是完整的扇形，而是以扇形兩側(cè)極點向上延伸，上方形成類似長方形結(jié)構(gòu)（圖8b）。

圖8 Poly-YOLO網(wǎng)絡(luò)分割過程

2.4 3種網(wǎng)絡(luò)模型分割性能

從表1可以看出，使用Poly-YOLO算法得到的TPR最高，HD最小，而使用馬爾可夫隨機場分割得到的DSC和IOU均最大。與Poly-YOLO和Ostu閾值分割相比，馬爾可夫隨機場在DSC的表現(xiàn)上分別提高了5%和3%，在IOU上分別提高了10%和2%。與Ostu閾值分割法和馬爾可夫隨機場相比，Poly-YOLO算法的TPR分別提高了11%和15%，HD分別降低了0.29和0.18。雖然對于扇形區(qū)域上方兩頂點的擬合程度較差，但Poly-YOLO算法仍然較為成功地完成了對于扇形區(qū)域的分割任務(wù)，在保證完整分割扇形區(qū)域的前提下，對周圍的非影像學(xué)信息進(jìn)行了剔除。圖9顯示了Poly-YOLO算法在各評價指標(biāo)的箱式圖中四分位數(shù)間距區(qū)間較窄，代表算法在面對不同結(jié)構(gòu)的肝囊型包蟲病超聲圖像時能表現(xiàn)出較為穩(wěn)定的性能，即具有較強的魯棒性，值得注意的是，在Ostu算法所得的箱式圖中，最下緣接近0，提示Ostu算法在處理數(shù)據(jù)集的某些圖像中，無法對扇形影像區(qū)域進(jìn)行較好地分割。

圖9 不同評價指標(biāo)的箱式圖

表1 不同分割算法性能比較

3 討論

本研究旨在測試Ostu閾值分割、馬爾可夫隨機場分割和基于深度學(xué)習(xí)的Poly-YOLO網(wǎng)絡(luò)模型分割3種方法在肝囊型包蟲病超聲圖像影像區(qū)域的分割性能。

Bindhu等[23]研究指出Ostu算法僅考慮灰度值，無法有效去除噪聲干擾，灰度直方圖易產(chǎn)生波峰波谷不明顯的問題，從而導(dǎo)致閾值計算錯誤，扇形區(qū)域難以被完整地分割，最終導(dǎo)致其在本研究中各項性能表現(xiàn)較差。如圖6d所示，在最大連通區(qū)域?qū)ふ业降难谀D中，Ostu閾值分割算法在處理邊緣模糊的區(qū)域時產(chǎn)生了過分割的現(xiàn)象，沒有將影像區(qū)域分割為完整的扇形，這種情況同樣出現(xiàn)在馬爾可夫隨機場分割算法中（圖7c），這可能是由于肝囊型包蟲病超聲圖像的特殊性，即病灶內(nèi)部為暗性液區(qū)[24]，導(dǎo)致算法無法對邊緣進(jìn)行有效的擬合，馬爾可夫隨機場應(yīng)用于圖像分割任務(wù)時，實質(zhì)上是將圖像分割問題轉(zhuǎn)化為對圖像中每一個像素進(jìn)行標(biāo)記分類的過程，通過求最大后驗概率來尋找每個像素的最優(yōu)標(biāo)簽，目標(biāo)點的標(biāo)記值僅與其鄰域內(nèi)的像素點標(biāo)記值有關(guān)，與鄰域外像素點無關(guān)，當(dāng)每個標(biāo)簽都得到了最合適的標(biāo)簽后，即代表分割完成[25]。

Fan等[26]研究指出馬爾可夫模型在MRI圖像分割任務(wù)上出現(xiàn)了由于灰度不均而導(dǎo)致模糊邊緣的過分割問題，扇形影像區(qū)域的不完整分割可能會對后續(xù)研究中關(guān)于肝囊型包蟲病亞型分類結(jié)果造成一定的影響。本研究發(fā)現(xiàn)，基于深度學(xué)習(xí)的Poly-YOLO算法可較好地處理這一情況，通過對大量標(biāo)注圖像的學(xué)習(xí)，該算法較為完整地實現(xiàn)了肝囊型包蟲病超聲圖像扇形影像區(qū)域的分割，這也使得其在TPR和HD指標(biāo)中具有較好的表現(xiàn)。但有研究指出，Poly-YOLO分割算法在處理非強凸圖像時，當(dāng)2個頂點同屬于一個極坐標(biāo)單元時，會取距離較大的頂點[14]，這可能是導(dǎo)致Poly-YOLO算法在預(yù)測扇形區(qū)域無法較好擬合扇形的2個頂點的原因，使其在DSC與IOU中表現(xiàn)不佳。

本研究仍存在一定的不足：馬爾可夫隨機場分割算法在圖像分割領(lǐng)域仍具有巨大的應(yīng)用價值，其利用圖像像素之間的關(guān)系，可較好地將屬于同一類型的像素精準(zhǔn)分類[27]。在本研究中，由于肝囊型包蟲病超聲圖像的特殊性，部分病灶與背景像素值相仿，導(dǎo)致算法將本該分為病灶的區(qū)域錯誤劃分為背景，致使其性能表現(xiàn)較差。未來研究可采用圖像增強相關(guān)算法，將目標(biāo)與背景區(qū)域?qū)Ρ榷仍鰪姡蕴岣唏R爾可夫隨機場對圖像分割的性能。在基于深度學(xué)習(xí)的Poly-YOLO算法中，由于網(wǎng)絡(luò)自身標(biāo)簽創(chuàng)建問題，導(dǎo)致模型對于扇形區(qū)域擬合并不完整，下一步的研究方向為對網(wǎng)絡(luò)模型及其參數(shù)進(jìn)行適當(dāng)?shù)恼{(diào)整，以便訓(xùn)練出的模型可以較好地擬合扇形區(qū)域，實現(xiàn)更為精確的分割效果。

4 結(jié)論

本研究結(jié)果表明，基于SSR的Ostu閾值分割方法、馬爾可夫隨機場圖像分割和基于深度學(xué)習(xí)的Poly-YOLO網(wǎng)絡(luò)模型均能較好地分割出扇形影像區(qū)域，其中Poly-YOLO算法相比較其他2種算法具有較強的魯棒性，適用于不同儀器報告的超聲圖像，且分割效果較好，TPR可達(dá)0.88。本研究為后續(xù)基于超聲影像的肝囊型包蟲病自動分型研究奠定了一定的理論基礎(chǔ)。