掃頻光學相干層析角膜圖像輪廓自動提取算法*

汪毅 劉珊珊 張瑋茜 蔡懷宇 陳曉冬

(天津大學精密儀器與光電子工程學院,光電信息技術教育部重點實驗室,天津 300072)

在掃頻光學相干層析系統中,遠心掃描模式造成角膜圖像中存在偽影、部分結構缺失及低信噪比區域,影響了角膜輪廓提取的精度.針對該問題,本文提出了一種針對低質量角膜圖像的輪廓自動提取算法.該算法首先依據圖像標準差分布將圖像劃分為高、低信噪比區域; 針對高信噪比區域,通過峰值點定位法獲取角膜輪廓; 針對低信噪比區域,通過連續幀圖像間配準疊加實現圖像增強,為低信噪比區域提供參考輪廓點,再通過權衡參考輪廓點與局部直線擬合結果的優劣,實現角膜輪廓定位; 最后,通過全局多項式擬合實現對全區域的角膜整體輪廓信息.對光學眼模型進行實驗,結果表明,與已有算法相比,本文算法對角膜輪廓的提取精度平均提高了4.9%.

1 引 言

掃頻光學相干層析成像(SS-OCT)是一種基于低相干干涉測量的醫學成像技術[1-3],具有非侵入、分辨率高、成像速度快的優勢,近年來,已成為臨床眼科檢查角膜的一種標準手段[4,5].SS-OCT系統獲取的角膜輪廓包含角膜曲率及角膜厚度等信息[6],這些信息是角膜異常定量分析及眼科臨床手術的重要依據[7-9].由于角膜圖像中常含有低信噪比區域與中心飽和偽影,同時存在輪廓信息缺失及模糊現象[10],而完整精準的角膜輪廓能夠提供更加可靠的病理信息,因此,實現對低質量的角膜圖像輪廓的精確提取十分必要.

傳統的角膜輪廓提取方法采用手動追蹤,但該過程冗長耗時且易出現人為誤差.為克服上述困難,多種自動提取算法被提出.Li等[11,12]提出結合快速活動輪廓和二階多項式擬合的方法,但該方法未考慮中央偽影且易受噪聲干擾; Shu和Sun[13]提出基于邊緣檢測和隨機抽樣一致的方法,雖然排除了中央偽影的影響,但未考慮角膜弱邊緣的存在;LaRocca等[14]提出魯棒性更強的基于圖論和動態規劃[15]的方法,雖然能夠實現對同時包含偽影及低信噪比區域的角膜輪廓的提取,但其對低信噪比區域的定位采用Li等[11,12]提出的二階多項式近似的假設,定位結果并不精確; Santos等[16]通過深度學習的方法實現對角膜的分割,但是該方法需要大量臨床數據及專業處理平臺.

本文提出一種針對SS-OCT系統獲取的低質量角膜圖像的輪廓自動提取算法,先根據圖像標準差將角膜圖像劃分為高、低信噪比區域并分別進行輪廓提取,再通過全局多項式擬合實現區域間的輪廓合并與重整.

2 算法原理

SS-OCT成像系統獲取的角膜圖像因其遠心掃描模式造成圖像中央區域信噪比較高,輪廓線條明顯,但含有偽影區域; 角膜兩側信噪比較低,角膜輪廓較為模糊,無法單純利用圖像梯度信息進行提取,同時SS-OCT系統的振鏡組高速掃描模式也造成角膜圖像輪廓及對比度并非連續均勻的現象.因此,針對包含偽影及部分結構缺失的高信噪比區域,可根據像素的灰度分布來定位峰值點作為角膜輪廓; 針對輪廓信息基本缺失的低信噪比區域,可通過幀間疊加對角膜圖像進行弱邊緣增強從而提供部分低信噪比區域的輪廓信息,同時利用該信息引導角膜輪廓局部直線擬合結果,使其更加貼近角膜真實輪廓; 針對角膜輪廓非連續均勻的問題,可通過全局多項式擬合實現對高、低信噪比區域角膜輪廓的整合及提取.角膜圖像輪廓算法流程如圖1所示.

圖1 角膜輪廓提取算法流程圖Fig.1.Flow chart of corneal contour extraction algorithm.

首先根據圖像標準差分布將圖像劃分為高、低信噪比區域; 之后對角膜圖像高信噪比區域進行輪廓提取; 接著利用連續多幀圖像疊加進行圖像增強,并通過最大類間方差法(OTSU)[17]及中值濾波算法實現對圖像背景噪聲的濾除; 之后通過權衡局部直線擬合結果及參考輪廓點實現對角膜低信噪比區域的輪廓提取; 最后,對高、低信噪比區域的所有輪廓點進行全局多項式擬合,實現對角膜輪廓的自動提取.

2.1 角膜高、低信噪比區域劃分

為實現對角膜圖像的分區處理,首先需要將其劃分為高、低信噪比區域.假設角膜圖像大小為M×N,M為A-scan的采樣點數,N為A掃描數,每個像素點的灰度值為I(i,j),其中i為像素點的行坐標,j為像素點的列坐標.利用(1)式計算各列像素間的標準差,為便于觀察,對(1)式的結果進行平滑濾波,處理結果如圖2(a)所示,其中橫軸信息為像素列坐標,縱軸信息為各列像素的標準差平滑結果.計算圖像各列像素標準差的均值meanall,并以meanall為閾值對圖像進行高、低信噪比區域的分割.在圖2(a)中搜索縱坐標大于meanall的區域,最終定位高信噪比區域為Region II,低信噪比區域為Region I,Region III.

2.2 高信噪比區域的輪廓提取

由于高信噪比區域輪廓線條明顯,其灰度明顯高于背景信號,一般表現為單列像素中的峰值像素點.因此對于高信噪比區域的輪廓,可根據其單列像素的灰度分布進行峰值點提取從而獲取角膜輪廓.由于角膜輪廓部分結構缺失及偽影區域的存在,導致單列像素峰值點提取過程中存在以下幾種情況：

圖2 角膜圖像高、低信噪比區域劃分 (a) 列間標準差平滑結果; (b) 角膜圖像高、低信噪比區域劃分結果Fig.2.Division between high and low SNR (signal-to-noise ratio) regions of corneal image：(a) Smoothing result of standard deviation between columns; (b) division results of high and low SNR regions of corneal image.

1)角膜上、下表面均不含偽影且不存在結構缺失;

2)角膜僅有一表面含有偽影;

3)角膜上、下表面均含有偽影;

4)角膜某一表面或上、下表面邊緣結構缺失.

為解決不同情況下對輪廓點的定位,首先設定圖像灰度均值threshold為判斷像素點是否為輪廓點的閾值,并依據人眼常見病變如高眼壓癥[18]、遠視眼[19]、近視眼[20,21]、青光眼[22]等以及人眼正常生理性變化[23],設定角膜上下表面輪廓點間距范圍為CT_range,偽影區域范圍為Artifact_range.

針對情況1),可直接搜索單列像素中大于threshold的峰值點及與峰值點間距大于CT_range的次峰值點作為角膜輪廓點.針對情況2),在搜索大于threshold的峰值點時,若同時存在多個像素點灰度相同或極為接近,且集中在Artifact_range范圍內,說明該像素點集處于偽影區域,則選取像素點集的中值點(針對灰度相同)或峰值點(針對灰度接近)作為輪廓點,其余點認定為偽影點,后續不予處理,并繼續在該列殘余像素中搜索與已定位的輪廓點間距大于CT_range的次峰值點作為另一表面輪廓點.針對情況3),在搜索大于threshold的峰值點時,若同時存在多個像素點灰度相同,且像素點中存在相鄰點對的間距超過CT_range時,則說明上、下表面輪廓點均存在偽影且偽影像素點與峰值點混疊,因此以相鄰點對為界將像素點集分為兩部分,并選取各部分中值點作為輪廓點,以保證偽影區域的輪廓點盡可能接近真實輪廓.針對情況4),當角膜結構缺失時,說明該列像素中不存在大于threshold的峰值點或存在峰值點但不存在與其間距超過CT_range的次峰值點,針對這種情況,可依據后續的局部多項式擬合算法進行缺失結構的填充.

根據上述方法對角膜進行初步輪廓提取的結果如圖3(a)所示,其中紅色點表示上表面輪廓點,綠色點表示下表面輪廓點.初步提取結果存在一定誤差,其主要來源于兩部分：上下表面灰度分布不均造成輪廓點定位結果與角膜實際結構不符; 孤立噪點引起的定位錯誤.

針對第一項誤差,可依據輪廓點定位結果的位置信息即坐標值大小將其劃分為上、下表面; 針對第二項誤差,可根據角膜結構具有局部穩定性[23]的特征,計算輪廓點中的所有相鄰像素的坐標差值,取其平均值為diff_mean,根據多次實驗結果,當某一輪廓點與其相鄰像素的坐標差值超過3diff_mean即可認為該點為誤差點.

經由上述算法的處理結果如圖3(b)所示,可以看出通過上述算法,高信噪比區域的角膜輪廓得到精確提取,偽影部分及缺損結構經由其相鄰像素的補充也得到了相對準確的定位.

2.3 圖像增強

單幀角膜圖像低信噪比區域與背景信號基本融為一體,信噪比極低,無法直接利用灰度分布來獲取角膜輪廓.由于SS-OCT系統成像速度極快,連續幀角膜圖像除少量剛性位移外,基本可以認為角膜結構不變; 而角膜圖像中存在包含散斑噪聲[24]在內的多種隨機噪聲,因此可以利用對噪聲容忍度較高的相位相關[25]算法對連續幀圖像進行配準疊加,幀間疊加過程中利用噪聲的隨機性對噪聲進行抑制同時實現對低信噪比區域角膜輪廓的增強.

經實驗驗證,對連續5幀角膜圖像進行配準疊加即可在保持角膜結構不變的前提下實現角膜圖像輪廓的增強.由于角膜圖像直接疊加會使得偽影區域大面積擴增并出現過飽和現象,造成偽影與角膜輪廓混疊以至于無法區分.因此首先利用2.1節分別對單幀圖像高、低信噪比區域進行劃分,接著利用2.2節算法分別對單幀圖像進行高信噪比區域的輪廓定位,最后對定位結果進行配準疊加,得到待處理的角膜增強圖像,實驗結果表明,經配準疊加后的角膜圖像低信噪比區域的信噪比從原來的單幀圖像的信噪比更高.

圖3 高信噪比區域的輪廓提取結果 (a) 輪廓點初步提取結果; (b) 輪廓點精確提取結果Fig.3.Contour extraction results of high SNR region：(a) Preliminary extraction result of contour points; (b) accurate extraction result of contour points.

2.4 圖像背景噪聲的濾除

由于角膜增強圖像中大部分像素為背景信息,因此可通過OTSU算法提取出角膜輪廓信息,提取結果如圖4(a)所示.由于上述濾波結果中仍存在孤立噪聲點,考慮到需要保持圖像的縱向輪廓信息,而角膜內部的細節信息在輪廓提取時并不重要.因此對圖4(a)結果進行橫向中值濾波,以實現對角膜圖像噪聲的抑制,濾波結果如圖4(b)所示.

2.5 低信噪比區域的輪廓提取

以高、低信噪比區域的邊界為起始搜索位置進行輪廓點的搜索,以待定位點的前向多個已定位輪廓點為基礎進行局部直線擬合,并以此預測該待定位點的坐標,同時權衡2.3節及2.4節提供的部分低信噪比區域參考輪廓結構信息與該預測值的優劣,確定最佳匹配輪廓點.設定最佳匹配輪廓點坐標為res,各像素點前向搜索范圍為range,以預測值為原點的角膜輪廓上下搜索范圍為2tolerance,經實驗驗證,本文選取的range=CT_range/4,tolerance=range/5.具體實現方法如下.

1)對第j列前向range列內的已定位輪廓點進行直線擬合,得到預測直線polyfit,由polyfit得到第j列的輪廓點預測值pred(j).

2)對低信噪比區域表面輪廓點進行建模(如圖5(a)所示),在第j列角膜各表面tolerance范圍內搜索非零像素點即參考輪廓點nonzero,若存在nonzero,則針對角膜上、下表面,分別進行以下判斷：

(a)對于上表面,若nonzero點坐標大于等于pred(j),即針對參考輪廓點,預測值更接近上表面,則pred(j)為最佳匹配輪廓點(見第11,10列); 若nonzero點坐標小于pred(j),則nonzero為最佳匹配輪廓點(見第8列);

圖4 角膜整體輪廓提取過程 (a) OTSU算法處理結果; (b) 中值濾波處理結果Fig.4.Extraction process of the overall cornea contour：(a) Processing result by OTSU algorithm; (b) processing result by median filtering.

圖5 低信噪比區域角膜輪廓建模及提取結果 (a)低信噪比區域表面輪廓點建模結果; (b)低信噪比區域輪廓提取結果Fig.5.Modeling and extraction results of corneal contour in low SNR region：(a) Modeling results of surface contour points in low SNR region; (b) contour extraction result of low SNR region.

(b)對于下表面,若nonzero點坐標小于等于pred(j),即針對參考輪廓點,預測值更接近下表面,則pred(j)為最佳匹配輪廓點(見第8,10列); 若nonzero點坐標大于pred(j),則nonzero為最佳匹配輪廓點(見第11列).

3)若搜索范圍內不存在nonzero,則obs(j)為最佳匹配輪廓點(見第9列).

通過上述方法實現對低信噪比區域的輪廓提取,最終得到的角膜邊緣如圖5(b)所示.

2.6 角膜完整輪廓提取

經過2.1-2.4節提取的高、低信噪比區域的輪廓點(圖6(a))并非平滑連續的角膜結構,因此需要對提取的輪廓點進行全局多項式擬合,通過整合高、低信噪比輪廓點的位置信息,以得到完整的角膜輪廓.擬合結果如圖6(b)所示.

圖6 角膜完整輪廓的提取過程 (a)角膜上下表面輪廓點提取結果; (b) 角膜輪廓擬合結果Fig.6.Extraction process of the complete cornea contour：(a) Extraction results of the contour points in the upper and lower cornea surfaces; (b) fitting result of the cornea contour.

3 實驗及結果

為驗證本文算法的正確性,首先搭建一套SSOCT成像系統對角膜進行成像.該系統選用中心波長為1060 nm,帶寬為30 nm,頻率為3 kHz,相干長度為69 mm的掃頻光源(AXSUN-1060),其在組織內的理論縱向分辨率為15 μm.樣品選用光學眼模型,其角膜折射率為1.376,角膜厚度為550 μm,符合正常成人角膜結構參數.對該樣品連續采集5幀圖像作為一組,共計采集10組圖像,采集的每幅圖像均包含820條A-scan,每條A-scan包含350個像素點,根據系統分辨率計算出角膜實際厚度對應的像素點數為Pstd.本文采用Intel Core i5-3337U處理器,4G內存,以Windows x64系統下的MATLAB 2017b為平臺進行實驗.

在現有算法中,只有文獻[14]能夠同時實現對含偽影及低信噪比區域的角膜圖像輪廓的提取,但其對低信噪比區域的處理采用傳統的二項式擬合[14],所獲取的弱邊緣輪廓精度較低.而本文算法利用了低信噪比區域的真實信息引導輪廓結構的提取,可增強提取精度與可靠性.

為驗證本文算法能否實現對低質量角膜圖像進行精確的輪廓提取,分別利用文獻[14]中的算法與本文算法對獲取的10組角膜圖像進行處理,兩種算法在實驗平臺上的平均數據處理時間分別為9.6 s和8.9 s.

處理的其中一組結果如圖7(a)所示,其他組實驗與該結果類似.圖中的紅色和玫紅色曲線為本文算法的提取結果,藍色和綠色曲線為文獻[14]提取結果,局部圖A為角膜上表面提取結果,局部圖B為角膜下表面提取結果,局部圖C為含偽影部分的高信噪比區域提取結果; 分別計算兩種算法獲取的角膜各位置的平均厚度,計算結果如圖7(b)所示,圖中藍色曲線表示本文算法獲取的角膜厚度值與Pstd間的偏差隨位置信息的變化情況,綠色曲線為文獻[14]獲取的角膜厚度值與Pstd間的偏差隨位置信息的變化情況,其中,紅色直線為角膜實際厚度Pstd; 分別計算兩種算法在高、低信噪比區域提取結果與Pstd間的平均誤差值與偏差率,結果如表1所示.

由圖7(a)的小圖A,B及圖7(b)可以看出,在低信噪比區域,本文算法提取的角膜輪廓更加貼近角膜上下表面.由圖7(a)的C可以看出,在高信噪比區域,本文算法與文獻[14]的提取結果相近,兩種算法相較于角膜的實際標準尺寸的平均測量偏差分別為7.2%和2.3%.本文算法獲取的角膜厚度更加接近角膜實際厚度,平均提取精度較文獻[14]的提高了4.9%.

圖7 兩種算法效果對比 (a)兩種算法輪廓提取結果; (b) 兩種算法角膜厚度平均計算結果Fig.7.Comparison of the effects of the two algorithms：(a) Results of the contour extraction of the two algorithms; (b) results of the corneal thickness calculated by the two algorithms.

表1 兩種算法對角膜輪廓平均提取精度對比Table 1.Comparison of the accuracy of two algorithmsfor contour extraction of high and low SNR regions.

4 總 結

本文提出了一種針對低質量的SS-OCT角膜圖像的輪廓自動提取算法,無需手動追蹤角膜輪廓,直接通過對高、低信噪比區域分別處理即可提取出低質量角膜圖像的角膜輪廓信息.該算法通過幀間疊加實現了圖像增強,繼而引入了角膜的實際輪廓信息作為權衡因子,增強了算法的魯棒性.與目前僅有的能夠處理既存在偽影又存在低信噪比區域圖像的角膜輪廓提取算法相比,本文算法對低信噪比區域角膜輪廓的提取精度更高,能夠為后續低質量角膜圖像輪廓的精確提取提供重要依據.