楊 昊 胡 曼 徐永利*
(1.北京化工大學 數(shù)理學院, 北京 100029; 2.北京兒童醫(yī)院 眼科, 北京 100045)
青光眼是世界第二大致盲性疾病,在2040年之前青光眼患者的數(shù)量會達到1.1億[1]。患者在青光眼的初期通常沒有明顯的癥狀,直到出現(xiàn)視覺損失時才會發(fā)現(xiàn)自己患病[2]。如果病人在青光眼初期就被發(fā)現(xiàn),采取及時的治療將能夠有效地阻止視力衰退。因此對青光眼的早期診斷十分重要。
視網膜神經纖維層(RNFL)缺損是診斷青光眼的重要特征。在臨床診斷中,光學相干斷層掃描(OCT)能夠直接測量RNFL厚度,被視為標準的檢查結果。但是OCT檢查價格昂貴,且我國多數(shù)中小型醫(yī)院和體檢中心并不具備相應設備。與OCT相比,眼底照檢查成本低,在許多醫(yī)院和體檢中心被廣泛使用。在眼底照上,RNFL于背景中呈現(xiàn)有條紋的模式,形成了特定的紋理。這些紋理會根據RNFL厚度的變化而改變[3]。因此,利用眼底照和OCT的多模態(tài)數(shù)據,設計基于眼底照來預測RNFL厚度的智能算法,不僅具有重要的臨床應用價值,也具有很高的可行性。
從20世紀80年代開始,許多基于眼底照提取RNFL特征的計算機輔助算法被研究者們相繼提出[3-9]。其中,Oliva等[8]提出了一種基于像素強度預測RNFL厚度的方法,并在9張正常眼和9張青光眼眼底照上進行實驗,Pearson相關系數(shù)達到了0.424。Odstrcilik等[9]在使用馬爾可夫隨機場、局部二進制等經典的圖像處理方法提取眼底照中的特征后,再通過支持向量回歸模型(SVR)預測RNFL厚度,并在28張正常眼和8張青光眼眼底照上進行評估;該方法在正常眼和青光眼上的Pearson相關系數(shù)分別為0.72和0.58,均方根誤差分別為26.06和25.09。然而以上文獻使用的傳統(tǒng)方法無法自動提取眼底照的RNFL信息,回歸預測值與真實值的偏差較大,而且沒有在大樣本數(shù)據集上進行驗證。近年來,深度學習在計算機視覺領域快速發(fā)展,并在醫(yī)學圖像領域取得了巨大進步[10-11]。姚紅革等[12]將深度學習方法應用在目標定位的回歸任務上,在公開數(shù)據集COCO 2014上達到先進水平。然而據我們所知,目前還沒有研究工作利用深度學習方法基于眼底照對視盤外圍局部區(qū)域的RNFL厚度進行預測。
本文設計了一種深度殘差回歸神經網絡,實現(xiàn)了自動提取眼底照中的特征并預測視網膜神經乳頭附近的RNFL厚度。利用本文設計的深度學習回歸算法,可以用較低的成本預測出眼底照中的RNFL厚度,量化評估眼底照中的青光眼性損傷,進而用于輔助青光眼的診斷。
圖1 OCT報告與眼底照配準示意圖Fig.1 Registration of OCT report and fundus image
本文所用數(shù)據來源于北京同仁醫(yī)院,均取自15 d內同時做過眼底照相、OCT檢查的正常人和青光眼患者的檢查結果,共計136只眼,來自于136個人,其中有107個正常眼和29個青光眼。
為了獲取OCT報告中RNFL厚度曲線在眼底照中所對應的位置,首先對彩色眼底照和OCT報告中的低分辨率眼底照進行圖像配準。本文使用PhotoShop軟件進行人工手動配準。在PhotoShop中將彩色眼底照置于底層,把這張眼底照對應的OCT報告中的低分辨率眼底照調整為半透明后置于頂層,通過調整低分辨率眼底照的位置和角度使兩者的血管重合,便完成配準。配準后即可得到OCT掃描路徑在彩色眼底照中的位置,記錄下掃描路徑所在圓周的圓心位置αi和半徑ri。
本文采用5折交叉驗證的方式進行實驗。在每折交叉驗證中,對于每張眼底照,在連續(xù)排列的5個patch中選擇1個劃入測試集,它的前一個劃入驗證集,剩下的3個劃入訓練集。劃分之后,訓練集共7 344個patch,驗證集共2 448個patch,測試集共2 448個patch,它們的組成如表1所示。在第g折交叉實驗中,第g組作為測試集,其余4組數(shù)據分別作為訓練集和驗證集(g=1,2,3,4,5)。訓練集用于訓練本文所提深度學習網絡以獲得預測模型,驗證集用于選擇網絡的超參數(shù),測試集不參與模型的訓練。最終,根據算法對5折交叉驗證實驗中測試集的預測來計算平均絕對誤差EMA(MAE)、決定系數(shù)和Pearson相關系數(shù)。
表1 訓練集、驗證集和測試集的組成Table 1 Components of train, validation and test dataset
基于眼底照中的正方形區(qū)域RGB圖像來預測該區(qū)域的局部RNFL厚度是一個回歸問題,因此我們設計了一種包含殘差模塊[14]的深度神經網絡來完成這一回歸預測過程。
殘差模塊最早出現(xiàn)于殘差分類網絡ResNet中,ResNet是ILSVRC2016比賽分類項目的冠軍,其跳躍連接結構提高了卷積層提取特征的能力并解決了深層神經網絡中的梯度消失問題[15]。但是相比于VGG16和Inception V3等網絡,經典的ResNet模型由于參數(shù)過多而使得模型復雜度較高。為了精簡模型參數(shù),本文利用全局平均池化層(GAP)來處理卷積層提取的特征。同時,將ResNet的損失函數(shù)由交叉熵損失替換為均方誤差損失,用于回歸預測任務。
以殘差神經網絡為主體結構,搭建了如圖2所示的深度殘差回歸神經網絡。網絡的輸入是大小為65×65的RGB圖像,輸入圖像經過卷積核大小為7、卷積核個數(shù)為64、步長為2的卷積層后,通過最大池化層減小維度,依次經過降采樣殘差模塊與殘差模塊。經過最后一個殘差模塊后連接全局平均池化層,后面加入全連接層,通過批標準化(BN)處理后以線性整流函數(shù)ReLU激活,再加入隨機失活率為0.7的DropOut層,最后使用只有1個輸出節(jié)點的全連接層作為回歸網絡的輸出。
圖2 深度殘差回歸神經網絡結構圖Fig.2 Architecture of the deep residual regression network
降采樣殘差模塊的詳細結構如圖3(a)所示,模塊的輸入經過一個ReLU激活層后有兩個分支,第一個分支依次經過3層卷積操作后,與第二個分支經過1個卷積操作的結果進行逐元素相加,作為降采樣殘差模塊的輸出。殘差模塊如圖3(b)所示,模塊的輸入經過一個ReLU激活層后有兩個分支,第一個分支依次經過3層卷積操作,之后與第二個分支進行逐元素相加作為殘差模塊的輸出。圖3中h1和h2表示卷積核個數(shù),s1和s2表示卷積操作的步長。每個階段中降采樣殘差模塊和殘差模塊的參數(shù)設置如表2所示。
圖3 降采樣殘差模塊和殘差模塊Fig.3 Sampling residual module and residual module
表2 深度殘差神經網絡中每個階段的參數(shù)設置Table 2 Parameters in each stage of the deep residual network
傳統(tǒng)的全連接層如圖4(a)所示,上一個卷積層輸出的特征圖譜伸展為一維向量后,經過全連接層與輸出節(jié)點相連。全連接層的參數(shù)量與這個向量的長度成正比,當特征圖譜維度過大時,全連接層的參數(shù)量激增,不僅增加了模型的體積也延長了計算時間。
圖4 全連接層與全局平均池化層Fig.4 Fully connected layer and global average pooling layer
損失函數(shù)選擇均方誤差損失EMS(MSE)來表示,計算公式為
(1)
式中,yi和y′i分別表示第i個樣本的真實值與算法預測值,mm。
選擇隨機梯度下降法作為優(yōu)化器,每次從全部樣本中隨機選擇一個樣本計算損失函數(shù)的梯度,沿著負梯度方向前進減小損失函數(shù)的值,通過不斷地更新權重使得損失函數(shù)值達到最小,此時的權重即為最優(yōu)。
選擇平均絕對誤差、決定系數(shù)R2和Pearson相關系數(shù)r作為算法的評估指標。這3種評估指標的表達式如下。
(2)
(3)
(4)
本文提出的深度學習回歸模型在以Tensorflow為底端的Keras框架下實現(xiàn),使用的設備是一臺裝有兩塊NVIDIA1080tiGPU的工作站。Tensorflow版本1.4.1,Keras版本2.0.1,Ubuntu系統(tǒng)版本16.04,內存64 G。
將原始圖像的像素值除以255縮放至[0,1]區(qū)間。本文采用遷移學習的方法,卷積層的初始參數(shù)通過對ResNet50在ImageNet數(shù)據集[16]上訓練獲得,并對全連接層的參數(shù)作隨機賦值;進而利用訓練集,通過最小化損失函數(shù)進一步調整網絡各層的參數(shù)。使用5折交叉驗證的方式對本文算法進行驗證,訓練過程中batchsize設置為64,迭代次數(shù)epoch設置為120,隨機梯度下降法的學習率設置為10-4。在訓練過程中選擇使驗證集損失最小的模型作為最終的模型,并驗證這一模型在測試集上的預測性能。
對本文所提算法在測試集上進行驗證,作為比較,也驗證了VGG16、InceptionV3、Xception和經典ResNet50(不含全局平均池化層)等網絡框架用于回歸預測時的預測性能,結果如表3所示。
表3 不同深度學習網絡在測試集上的評估指標Table 3 Performance of different deep learningnetworks on the test set
從表3可以看出,與VGG16、InceptionV3和Xception相比,經典ResNet50具有更小的MAE、更大的R2和r,但是模型體積增加較大。與經典ResNet50相比,本文提出的加入GAP的殘差網絡模型體積減小了19.6%,預測精度與其大致相當:MAE略有減小,R2略有降低,r略微提升。從綜合性能考慮,本文所提方法的預測效果最佳。
圖5展示了本文方法對RNFL厚度的預測值與OCT測量值的分布,可以看到兩者之間存在強相關性(r=0.878,P<0.001)。算法預測的RNFL厚度均值為(88.45±1.39) mm,OCT測量的RNFL厚度均值為(89.69±1.67) mm。進一步,分別在正常眼和青光眼上計算了算法預測值與OCT測量值的相關性,得到對于正常眼r=0.885,對于青光眼r=0.872。
圖5 算法預測值與OCT測量值的分布Fig.5 Distribution of prediction and OCT
為了分別觀察正常眼和青光眼的絕對誤差分布,圖6展示了在兩組不同數(shù)據上絕對誤差的密度分布曲線。可以看出,正常眼與青光眼的絕對誤差的密度分布曲線比較相似,誤差值大部分都分布在20 mm以下,其中正常眼在[0,20]區(qū)間的比例更高。正常眼的預測絕對誤差為(14.884±12.522) mm,青光眼的預測絕對誤差為(15.108±13.215) mm,青光眼的值略高,但差異不明顯。
圖6 絕對誤差在正常眼和青光眼上的分布Fig.6 Distribution of absolute error in normal and glaucomatous eyes
為了觀察絕對誤差在360°范圍內不同角度處的分布,將眼底照以配準后所得圓心為中心,從0°開始等分出10個扇形組,將每個扇形組的絕對誤差以箱線圖的形式展示在圖7中。在箱線圖中,矩形的水平上邊框和水平下邊框分別表示絕對誤差的3/4分位數(shù)和1/4分位數(shù),矩形內部的水平直線表示中位數(shù),矩形上方的“+”表示奇異值??梢钥闯?,多數(shù)扇形區(qū)域的絕對誤差的3/4分位數(shù)都分布在25 mm以下,每個扇形區(qū)域內都有少數(shù)奇異點大于30 mm。其中第2組和第9組的箱線圖整體分布的位置更高,且3/4分位數(shù)超過25 mm,遠高于其他組。這兩組分別對應36°~72°和288°~324°(視網膜神經乳頭左上側和左下側),這兩處的血管較粗且分布較為密集,影響了模型在附近區(qū)域提取RNFL特征的性能。然而,即使是不同的OCT設備,測量出的RNFL厚度值也會存在10~25 mm的差異[9],因此本文方法的誤差是可以接受的。
圖7 算法預測值的絕對誤差在不同角度的分布Fig.7 Distribution of the predicted absolute errors at different angles
進一步在整張眼底照上進行了算法預測值與OCT測量值的對比。圖8展示了在一張正常眼眼底照上本文算法預測的RNFL厚度值和OCT測量值的對比,可以看到正常眼的預測曲線與OCT測量曲線在各個位置都很接近,兩者均刻畫出了正常眼的典型雙峰結構。圖9展示了在一張青光眼眼底照上本文算法預測的RNFL厚度值與OCT測量值的對比,同樣可以看到預測曲線與OCT測量曲線在絕大多數(shù)位置都很接近,而且在70°~170°之間準確預測出了此處的RNFL厚度接近于0,表明此處存在RNFL缺損。
圖8 正常眼上算法預測值與OCT測量值的對比Fig.8 Comparison between the prediction and OCT on a normal eye
圖9 青光眼上算法預測值與OCT測量值的對比Fig.9 Comparison between the prediction and OCT on a glaucomatous eye
本文基于眼底照和OCT的多模態(tài)數(shù)據,提出了基于眼底照中的局部區(qū)域信息預測RNFL厚度的智能算法。我們設計了以殘差卷積結構提取圖像特征、以全局平均池化層與全連接層輸出回歸預測值的深度殘差回歸神經網絡。在真實數(shù)據集上的實驗結果驗證了本文算法具有較低的平均絕對誤差和較高的統(tǒng)計相關性。通過本文算法,可以基于成本較低的眼底照來預測RNFL厚度,也可以與其他青光眼類指標聯(lián)合來輔助青光眼的診斷。