馬嘯天,李書芳
(北京郵電大學信息與通信工程學院,先進信息網絡北京實驗室,網絡體系構建與融合北京市重點實驗室,北京 100876)
基于深度學習方法的核磁共振圖像(Magnetic Resonance Image,MRI)分割是圖像分割研究領域中的研究熱點,在醫(yī)學領域具有重要的應用價值。心臟是人體最重要的器官之一,是循環(huán)系統(tǒng)中的主要器官,為血液流動提供動力。心血管疾?。–ardiovascular Disease,CVD)具有高患病率、高致殘率和高死亡率的特點,已經成為人類最普遍的死亡原因[1]。心臟功能分析在臨床心臟病學中對患者管理、疾病診斷、風險評估和治療決策發(fā)揮著重要作用??梢酝ㄟ^定量評估心室容積、卒中體積、射血分數和心肌質量等臨床參數,對心臟功能進行分析。這些參數的計算前提是對左心室、右心室和心肌輪廓的精確分割。目前臨床上評價心臟功能的影像技術中,心臟磁共振成像(Cardiac Magnetic Resonance Imaging,CMRI)技術良好的時間和空間分辨率使其具備同時顯示心臟結構和功能的能力,加之其不存在輻射損害,這種集形態(tài)、功能及細胞生物學檢查為一體的無創(chuàng)性檢查已發(fā)展為心臟病診斷和鑒別診斷的理想方法,被認為是判斷心臟結構和功能的“金標準”[2]。因此心臟MRI 分割已成為一個醫(yī)學影像分析的熱點問題。
通常,放射科醫(yī)生通過手動繪制包含分割目標結構的輪廓線,從周圍的組織器官中劃定目標區(qū)域。由于醫(yī)學影像本身的復雜性以及對結果極高精確度的要求,目前醫(yī)學影像的分析工作主要由有經驗的醫(yī)生完成,自動化的分割與診斷由于不能滿足臨床需要,只能作為輔助補充。然而放射科醫(yī)生人工分析所需的工作量大、耗費時間長,且受不同的放射科醫(yī)生各自主觀的經驗、所處環(huán)境、工作狀態(tài)等影響,其結果也因人而異,對同一區(qū)域的勾畫結果不能百分之百復現。隨著醫(yī)療數字化的推進,利用各種設備采集到的醫(yī)學圖像數量也快速增長,這些更增加了醫(yī)生的工作量。在現有的醫(yī)學影像分割技術中,深度學習方法體現出了明顯的優(yōu)勢。
將深度學習技術用于計算機輔助診斷中不僅可以提高診斷精度和效率,還能有效緩解醫(yī)療資源緊張,為醫(yī)生減輕壓力,降低對專家經驗的依賴,提升基層診療水平。然而在醫(yī)學影像分割任務中,也面臨著諸多挑戰(zhàn):由于掃描成本、患者隱私以及高人工標注成本等問題,使得用于訓練深度學習模型的數據規(guī)模十分有限;醫(yī)學圖像的輸入維度高,通常為體積數據,數據量大;需要分割的目標體積較小,對精度的要求更高;不同的分割目標差異性較大,缺少可遷移到醫(yī)療影像分割任務上的通用預訓練模型等。心臟MRI 自身的特點,例如左心室腔內亮度不均勻,與其他器官亮度信號相似,心臟病引起心室心肌形狀畸形,以及心尖和基底切片圖像的復雜性等等,也給心臟MRI 的分割任務帶來諸多困難。
基于上述問題,本文提供了基于深度學習的心臟MRI 分割方法。通過感興趣區(qū)域檢測及深度學習網絡技術,對心臟MRI 圖像與醫(yī)生手動分割的真值標注數據進行訓練,實現對心臟MRI 圖像左右心室和心肌的準確、自動分割。使用了多種數據增強手段進行實驗數據集的人為擴充。Jaccard 系數,Dice重疊系數和豪斯多夫距離(Hausdorff Distance,HD)作為分割效果的客觀評估指標[3]。與現有的幾種方法相比,我們的方法取得了優(yōu)秀的結果。
過去幾年來,心臟MRI 的左心室、右心室的心內膜和心外膜分割是該領域的研究重點,為了推進先進的醫(yī)學方法也已經組織了一些半自動/全自動的心臟分割挑戰(zhàn),這些挑戰(zhàn)通常提供由專家標注的真實的數據集,并提供一組評估指標來對各種方法進行評價,如2009年和2011年的左心室分割挑戰(zhàn)[4,5],2012年的右心室分割挑戰(zhàn)[6]和2015年的Kaggle 第二屆年度數據科學碗挑戰(zhàn)賽,2017年MICCAI 的自動心臟診斷挑戰(zhàn)ACDC[7]等。
本文我們使用ACDC 數據集。數據集由150位不同患者的心臟MRI 組成,包括100個訓練樣本和50個測試樣本,每個訓練樣本包含舒張末期(End of Diastole,ED)和收縮末期(End of Systole,ES)的左心室、右心室以及心肌(心外膜輪廓)的專家手動分割標注結果。每位患者的MRI 數據包含28-40 幀從左心室底部到頂部的整個心動周期的一系列短軸圖像切片。每個切片的空間分辨率平均為235-263個體素。
Petitjean 和Dacher(2011)[8]對先前的左心室分割方法進行了全面的闡述,包括多級Otsu 閾值、可變形模型和水平集、圖形切割、基于知識的方法(如主動和外觀形狀模型)和基于Atlas 的方法。這些方法在小基準左心室數據集上取得了一定的成功,但它們的魯棒性差、準確性低,且泛化能力有限。2015年Kaggle 挑戰(zhàn)賽的結果表明,性能最高的方法依賴于深度學習技術,尤其是全卷積網絡[9]和U-Net[10]。Avendi 等人(2015)[11]將機器學習的最新進展(例如堆疊自動編碼器和卷積神經網絡)與可變形模型相結合,使左心室心內膜分割的結果達到了新的水平。2016年,P.V.Tran[12]首次使用全卷積神經網絡同時分割出左心室和右心室的內外膜。2018年,Nasresfahani 等[13]首先提取出感興趣區(qū)域(Region of Interest,ROI),在此基礎上采用全卷積神經網絡,提高了分割精度,同時提高了模型的魯棒性。
本文提出的心臟MRI 分割方法由數據預處理與分割兩部分組成。
預處理步驟包括對心臟MRI 數據的感興趣區(qū)域檢測、數據標準化、數據增強。
3.1.1ROI 檢測
ROI 檢測的步驟包括使用標準差濾波器進行濾波,找到心動周期圖像序列中像素強度隨時間變化強烈的區(qū)域,然后使用Canny 邊緣檢測[14]和霍夫變換技術[15]找到包含心臟的感興趣區(qū)域。
Canny 邊緣檢測算法由多個步驟構成:
(1)圖像平滑。采用高斯平滑濾波來去除圖像中的噪聲。在圖像邊緣檢測過程中,圖像的邊緣和噪聲難以區(qū)分,僅憑邊緣檢測算法不能完全消除噪聲對邊緣檢測過程和結果的影響,因此需要對原始圖像進行預處理。圖像預處理中的常見濾波方法有均值濾波、中值濾波和高斯濾波。相較于均值濾波和中值濾波,高斯濾波在對圖像進行平滑濾波時能很好地保留圖像中的灰度分布情況。令f(x,y)為輸入數據,fs(x,y) 為經過高斯卷積平滑后的圖像,G(x,y)表示高斯函數,表達式如下:
(2)計算圖像梯度強度和方向。Canny 算法的基本思想是尋找一幅圖像中灰度強度變化最強的位置,即梯度方向。平滑后的圖像中的每個像素點的梯度由Sobel 算子來計算。首先運用如下的卷積陣列Sx和Sy分別求得原始圖像A沿水平(x)和垂直(y)方向的梯度Gx和Gy:
然后利用下列公式求得每一個像素點的梯度幅值:
在變化劇烈的地方(邊界處)將獲得較大的梯度度量值G,然而這些邊界通常非常粗,難以標定邊界的真正位置,為了標定邊界,還需要梯度的方向信息:
(3)應用非極大值抑制(Non-Maximum Suppression)技術來消除邊緣誤檢。對每個像素點將其梯度方向近似為(0°,45°,90°,135°,180°,225°,270°,315°)中的一個,比較該像素點和其正負梯度方向的兩個像素點的梯度強度,如果該像素點梯度強度最大則保留,將其視為邊緣,否則抑制。經過這樣的處理后會保留一條邊界處最亮的一條細線,圖像的邊緣明顯變細。但是這些細線部分可能是噪聲或其他原因造成的局部極大值。
(4)應用雙閾值來決定可能的邊緣。Canny 算法中應用了雙閾值的技術,即設定一個閾值上界T1和閾值下屆T2,如果圖像中的像素點梯度值超過T1,稱為強邊緣,小于T1大于T2的稱為弱邊緣,小于T2的不是邊緣。通過設置足夠高的T1,強邊緣像素點的梯度足夠大(變化足夠劇烈),強邊緣必然是邊緣點。為了判斷弱邊緣是邊緣還是噪聲,當弱邊緣的周圍8鄰域有強邊緣像素存在時,就將該弱邊緣點保留為邊緣像素,以此來實現對強邊緣的補充,若沒有,則去除。Canny 建議高閾值與低閾值T1:T2之比為2:1。
對于霍夫變換,霍夫變換的原理是運用兩個坐標空間之間的變換,利用點和線的對偶性,將一個空間中具有相同形狀的直線或曲線映射到另一個坐標空間中的一個點上形成峰值,從而把檢測任意形狀的問題轉化為統(tǒng)計峰值問題。
以圓形霍夫變換(Circular Hough Transform,CHT)為例,霍夫變換首先提取圖像的邊緣,然后利用圓的方程將邊緣圖像中的每個點映射到另一個三維空間中去,這個三維空間稱為霍夫空間?;舴蚩臻g的前兩個維度表示待檢測的圓的圓心坐標,第三個維度表示待檢測的圓的半徑。邊緣圖像中的每個點在霍夫空間中都被映射為一個圓錐面,同一個圓上的點在霍夫空間中的映射會相交于一點,該點對應的坐標就是待檢測的圓的參數。
二維空間中的圓的方程用下式表示:
圓上任意一點的(xi,yi)被映射到三維參數空間a-b-r中的一個曲面上,該曲面方程用下式表示:
將圓上的每個點都進行以上的映射過程。假如每一個點對應的圓錐面都相交于點()a0,b0,r0上,則所有的點都在這三個參數所確定的圓上,這個圓就是要檢測的圓。
本文在Canny 邊緣檢測的基礎上,使用霍夫變換進行圓環(huán)的檢測來進行ROI 檢測,找到左心室區(qū)域的中心。
3.1.2數據標準化
由于采集設備和條件的不同,不同患者之間的心臟MRI 圖像可能由于明暗、對比度等原因,出現尺度不一的問題。而神經網絡學習的本質上就是為了學習輸入數據的分布,如果訓練數據和測試數據的分布差異過大,會大大降低模型的泛化能力。為了解決這個問題,通常需要對輸入數據進行標準化處理。
將采集到的心臟短軸MRI 數據轉化為像素灰度均值為0,方差為1 的標準化數據;使用Z-score 標準化方法[16]對圖像數據的像素強度進行標準化,
其中X為切片圖像像素,μ為圖像像素均值,σ為圖像像素標準差。經過標準化處理后的圖像像素值分布符合標準正態(tài)分布,均值為0,方差為1。
3.1.3數據增強
由于掃描成本、患者隱私、高質量的圖像標注獲取困難等原因,大量基于深度學習的醫(yī)學圖像研究面臨可用數據集少、訓練數據量不足的問題。小樣本數據集往往不足以支撐深度神經網絡訓練,常出現模型在訓練樣本上表現很好但在測試數據集上泛化效果不佳的過擬合現象。為了解決訓練數據量不足的問題,本文采用數據增強技術來人為地增加數據集,以防止模型過擬合。常見的數據增強方法有縮放、翻轉、旋轉、裁剪、平移等幾何變換方法和添加噪聲、模糊、顏色變換、擦除、填充等顏色變換方法。本文采用的數據增強方法包括-5~5 度的隨機角度旋轉,-5~5 mm 的隨機平移,-0.2~0.2 倍數的隨機縮放比例,添加零均值的標準差為0.01 的高斯噪聲,以及使用B 樣條插值進行彈性形變。
3.2.1網絡結構
本文提出的分割網絡結構基于Jégou 等人2017年提出的語義分割架構[17]進行改進,改進后的網絡結構如圖1 所示。網絡結構類似U-Net 的典型的語義分割結構,包括下采樣和上采樣兩個部分,由7個密集連接塊[18],3個過渡層上采樣,3個過渡層下采樣組成。在下采樣路徑中,密集連接塊的輸入與輸出通過特征通道進行拼接。在上采樣路徑中,通過特征圖元素級相加,與對應的下采樣路徑中的特征圖進行融合,其中改進的跳躍連接使用線性投影操作來匹配對應的特征通道數。
圖1 網絡結構
在輸入128×128×1 的圖像后,首先通過改進的Inception 模塊[19]獲得圖像的多尺度特征。其中原始的Inception 模塊由步長為2 的最大池化操作和卷積核大小分別為1×1、3×3、5×5 的卷積操作組成,改進后的Inception 模塊由3×3、5×5 和7×7 大小卷積核的卷積操作組成,去除了最大池化操作,并引入了更大的7×7 卷積核的卷積操作,擴大了模型的感受野,有利于大目標區(qū)域的檢測,并減少與目標相似的假陽性區(qū)域的誤判。
下采樣部分類似于DenseNet[18]的結構,在上采樣部分通過轉置卷積、密集連接和改進的跳躍連接[10],使下采樣的結果恢復到原始輸入的空間分辨率。改進的跳躍連接由批次歸一化(Batch Normalization,BN)、ReLU 激活函數、1×1 卷積操作和dropout 組成,使得下采樣和上采樣過程中不同維度的特征通道通過投影操作相匹配,實現了跨通道的特征信息交互。使用元素級相加進行拼接,與直接的特征通道拼接相比,減少了參數量和內存占用。
最后將上采樣得到的特征圖與1×1 卷積層和softmax 層相結合,使用Argmax 生成分割的最終結果圖。
3.2.2損失函數
在心臟MRI中,心臟相關結構對應的像素為前景像素,心室周圍組織器官對應的像素為背景像素。在整個心臟MRI 中,絕大部分的像素都屬于背景像素,因此網絡預測的結果可能會更偏向于背景類像素,如果不給前景像素設置足夠的權重,會導致網絡訓練傾向于給背景類像素更高的權重,造成大量像素的錯誤分類。在醫(yī)學圖像中,這種嚴重的類不平衡問題通常由損失函數解決,例如Dice 損失[20]和加權交叉熵損失[21]等。
交叉熵損失通過計算預測的輸出類和目標類之間的體素誤差概率來度量所有體素的累積誤差。假設y是真實的概率分布,是預測的概率分布,即神經網絡的輸出,則交叉熵函數的公式定義為:
其中yi表示樣本i的指示變量,如果樣本i的預測類別與該樣本的類別相同則值為1,否則為0;表示樣本i類別的預測概率。
如果有N個訓練樣本,則整個樣本的平均交叉熵為:
其中M表示類別的數量;yi,c表示樣本i的指示變量,如果樣本i的預測類別與該樣本的類別相同則值為1,否則值為0;表示觀測樣本i屬于類別c的預測概率。
考慮單獨使用交叉熵損失無法解決這種出現的類不平衡問題,本文設計加權方案用于加權交叉熵損失的計算。
從手工分割真實標注圖像生成權重圖,用于加權交叉熵損失公式在每個體素處計算的損失。加權交叉熵損失LW由下式計算:
其中,W=(w1,w2,…,wl)為可學習權重的集合,wl為與深度學習網絡第l層對應的權重矩陣,X為訓練樣本,ti為體素xi?X對應的目標類標簽,w(xi)為在每個體素xi處估計的權重,p(ti|xi;W)表示網絡最后輸出層的對體素xi的概率預測。
V-net[20]提出了圖像分割領域中另一個常用的損失函數Dice-loss。Dice 損失函數是根據dice 系數設計的。V-net中的損失函數為:
其中pi是第i個像素值預測的結果,gi是其真實值,1 代表目標,0 代表背景。Dloss的取值在0 到1 之間,越接近1表示預測的結果越接近真實值。
最終損失函數使用加權交叉熵和Dice 系數的組合,并使用L2正則化,形式如下式所示:
其中,λ、γ和η是單個損失的權重參數,||?||為L2范數,||W||2為L2 正則化項,其權重參數設置為η=5× 10?4。
本文網絡基于Python 3.6 和Tensorflow 1.12 實現。實驗均在服務器的單個Nvidia Geforce GTX 1080ti GPU卡上進行。網絡訓練使用ADAM優(yōu)化器,學習率設置為0.001,損失函數由式(14)確定。數據集使用MICCAI 在2017年提出的ACDC 數據集,共150例不同患者的心臟短軸電影序列MRI 圖像數據,其中包含有左心室、右心室以及心?。ㄐ耐饽ぽ喞┑膶<沂謩臃指顦俗⒌牟±?00例。網絡訓練批次由16幅預處理后的128×128 的二維MRI 圖像組成,設置訓練epoch 為250。每個epoch 后,在驗證集上對模型進行評估,確保最終選擇的用于測試集評估的最佳模型在驗證集上具有最佳的分割評估效果。
讀取原始心臟MRI數據,經過包含ROI檢測的預處理步驟后,得到結果如圖2 所示,其中(a)為原始心臟MRI數據輸入,(b)為經過ROI檢測后,以左心室為中心繪制了包含整個左心室的ROI灰色掩膜。
圖2 ROI檢測結果
經過圖像數據標準化和數據增強處理后,根據ROI 檢測得到的ROI 區(qū)域中心和區(qū)域半徑,將圖像數據裁剪為以ROI 中心即左心室為中心的128×128大小的塊,作為深度學習分割網絡的輸入。與將每張切片空間分辨率平均為235-263個體素的原始心臟MRI 切片數據直接作為輸入相比,使用相同的模型訓練占用的GPU 顯存大小從10GB 以上下降為不到6GB。
ROI 檢測分割結果如圖3 所示,圖3 中(a)為輸入圖像,(b)為經過ROI裁剪后的128×128大小的圖像。
圖3 ROI檢測分割結果
將經過預處理后的圖像數據輸入深度學習分割網絡訓練250個epoch 后,將訓練得到的模型應用到測試集上得到分割結果。以一例患者心臟ED 時的分割結果為例,如圖4 所示,圖中(a)為患者心臟ED 時的原始MRI,(b)為ED 時的專家手動分割標注,(c)為ED 時的網絡分割結果,(b)和(c)中白色為左心室區(qū)域,淺灰色為心肌區(qū)域,深灰色為右心室區(qū)域。可以觀察到本文所提出方法分割結果與專家手動標注結果非常接近,準確分割了左心室、右心室和心肌。
圖4 患者分割結果
本文采用Jaccard 系數,Dice 重疊系數和豪斯多夫距離來評估測試數據分割效果,它們在醫(yī)學圖像處理的評估上已經得到標準化。
Jaccard 系數又稱為Jaccard 相似系數,用于比較樣本集之間的相似性與差異性,Jaccard 系數值越大,樣本相似度越高。Jaccard 系數定義如下式,其中P和G分別為預測和真實標簽輪廓所包圍的像素集:
Dice 系數基本上衡量了預測的標簽與真實值的相似程度,Dice系數取值為0到1之間,越高的Dice系數表示分割的性能越好。Dice系數定義如下式:
豪斯多夫距離是描述兩個點集之間距離的一種度量,度量了兩個點集間的最大不匹配程度,高豪斯多夫距離表示兩個點集不緊密匹配。豪斯多夫距離定義如下式:
其中p和g分別為點集P和G中的點,d(p,g)為兩點間的歐幾里得距離,||?||為L2范數。
本文分割模型訓練使用的損失函數由式(14)確定,式中λ、γ和η均設置為1。與使用其他損失函數訓練作對比,對分割結果分別使用Dice 重疊系數、HD、Jaccard 系數進行評估,結果見表1、表2 和表3,其中CE 表示交叉熵損失函數,Dice 表示Dice 損失函數,SwCE 表示加權交叉熵損失函數,L2 表示L2 正則化,ED 和ES 分別表示心臟心動周期的舒張末期和收縮末期。表中數值為平均值,括號內為標準差。
表1 不同損失函數分割結果評價(Dice系數)
表2 不同損失函數分割結果評價(HD)
表3 不同損失函數分割結果評價(Jaccard系數)
將分割結果評估前三用粗體顯示,可以看到本文提出的加權交叉熵損失函數SwCE 有顯著的優(yōu)越性,在左心室、右心室和心肌上的分割結果評估均獲得領先,證明本文提出的加權方案可以更好的解決醫(yī)學圖像中的類不平衡問題。
使用由式(14)確定的SwCE、Dice、L2 正則化相結合的損失函數的分割結果,在豪斯多夫距離上取得了最優(yōu)的評價,表明分割的整體效果接近專家手動標注結果,在心臟的基底和心尖也沒有出現偏差過大的誤判,分割模型的魯棒性要比單獨使用Sw-CE 損失函數更好。
我們比較了本文提出的模型與2017年ACDC 挑戰(zhàn)前三名的模型,在50例測試集上的分割結果評估對比如表4和表5所示。結果表明本文提出的方法與挑戰(zhàn)前三名獲得了相近的結果,并且在HD 這一評估指標上結果優(yōu)秀,表明與專家手動標注結果更接近,誤判區(qū)更少,在心臟的基底和心尖切片處的分割效果更穩(wěn)定,模型的魯棒性更好。
表5 與其他模型對比(HD)
本文提出了一種基于深度神經網絡的心臟MRI自動分割方法。提出的心臟MRI 分割方法由預處理和分割兩部分組成。預處理步驟包含了使用Canny邊緣檢測和霍夫變換等方法的感興趣區(qū)域檢測、數據標準化、數據增強等,能夠減少模型訓練占用的內存大小以及非目標域信息帶來的干擾,改善訓練數據不足時分割模型的性能。預處理后,我們對檢測到的ROI 進行心臟結構的自動分割,提出的分割模型基于Jégou 等人2017年提出的語義分割架構進行改進,模型的訓練使用提出的加權交叉熵、Dice 系數和L2 正則化組合的損失函數。我們在ACDC 數據集上實驗了我們的方法,并與其他人作對比,實驗結果表明,我們的方法的分割結果接近專家手動分割結果,在客觀評價指標上獲得了優(yōu)秀的表現,且該方法具有更好的魯棒性。