基于圖搜索和梯度弱化的SDOCT圖像ILM層分割

摘 要： 頻域光學斷層相干掃描（SDOCT）圖像針對視網膜中的RNFL層的精確檢測，對于一些眼科疾病診斷有很好的參考價值，需要精確分割該層的上表面ILM邊界。給出一種3D圖搜索與梯度弱化方式相結合的分割ILM表面的方法，提出一種梯度弱化的方法，弱化了圖像中其他高對比度區域的梯度，降低其他層對ILM表面分割的干擾，從而在圖像中僅有在ILM表面附近較大梯度的前提下，利用經典三維圖搜索方法實現對ILM表面的精確分割。分割過程中采用多尺度的分割方法，降低了三維圖搜索算法本身時間復雜度和空間復雜度。采用斯坦福大學的25對眼睛的頻域光學相干斷層（SDOCT）圖像，共50組3D圖像體數據實驗，結合ITK?SNAP軟件的半自動分割結果進行實驗對比。實驗結果表明給出的方法對50組數據中的43組數據能夠實現ILM表面的快速精確分割。

關鍵詞： 3D圖搜索； 梯度弱化； 多尺度分割； ILM表面分割

中圖分類號： TN919?34； TP391 文獻標識碼： A 文章編號： 1004?373X（2015）04?0106?05

0 引 言

OCT圖像層分割在常見的眼科疾病中（青光眼、AMD（Age?Related Macular Degeneration，老年黃斑性變性））對于病變區域的診斷有很高的參考價值 [1]。圖1給出了一幀黃斑中央凹附近區域單幀的圖像分層示意圖。

近年來針對視網膜分層的方法主要包含如下幾類。文獻[2]提出基于區域增長的分割方法，該方法不能保證所求表面全局最優；Kajic等人給出了一種基于圖像紋理和形狀統計信息的邊緣特征提取分割算法，醫學圖像的紋理特征提取相對比較困難[3]； 文獻[4]采用結合形態學的邊緣提取方法，該方法針對部分邊緣不明顯的醫學圖像效果不甚好，沒有利用幀間信息；文獻[5]給出了基于水平集方法分割OCT圖像，其對應的符號距離函數的計算量較大。以上幾類分割方法均以OCT圖像的一個B?Scan（一幀切片圖像）為單位進行圖像分割，沒有充分考慮連續圖像幀與幀之間相對連續關系；目前基于圖論的分割方法可分成2D和3D兩類。2D的圖論算法的圖像分割單位針對單幀圖像。文獻[6]采用圖搜索的方法對OCT視網膜圖像進行層分割，利用Dijkstra算法搜索出最短路徑，逐步區域限制分割出各個視網膜層，其針對一些邊緣不明顯圖像也需要做部分分割后處理操作。基于3D圖論的方法在近年來得到廣泛的應用。3D圖搜索采用構圖的方法，給出了多表面同時分割的理念并引入3D子圖間限制，最終轉換成求解一個4D結構圖的最優解問題，目前Dufour等人在3D圖搜索基礎上又給出了先驗信息模型的概念，并結合先驗信息模型給出了軟限制構圖的思想，但依然需要層與層之間約束構圖同時分割求解，不能做到簡單，單獨分割某層[7]。如圖1所示。其中表面標記從上到下分別為內界膜（Internal Limiting Membrane，ILM），INL?OPL（Inner nuclear layer? Outer plexiform layer），OPL?ONL（Outer plexiform layer? Outer nuclear layer），IS?OS（photoreceptor innerand outer segments），OS?RPE（outer segments?retinal pigment epithelium）和BM（Bruch's membrane）。

直接使用3D圖搜索算法求解較慢，并且難以單獨分割某一層。本文在結合3D圖搜索方法基礎上，引入一種梯度弱化的方法，并將該方法應用到ILM表面的分割問題上。該方法優勢是，第一，能夠單獨分割ILM層。直接單獨分割ILM層可能會產生誤分割，常采用3層同時分割才可分出ILM層，當僅需分割ILM層時不甚合理；第二，利用梯度弱化后結合3D圖搜索分割ILM層會加快圖搜索的速度，原因是僅保持ILM表面附近較大梯度，其他區域梯度較小且相對一致，頂點權重w的計算方式導致構圖過程很多頂點權重為0，在將3D圖搜索轉化成求解最小割問題后，這些點會排除在s?t圖以外，大大減少s?t圖中的邊數，加快求解最小割速度。

1 方 法

本節給出梯度弱化方法的概念及算法實現。算法實現流程如圖2所示。

1.1 3D圖搜索算法模型

3D圖搜索算法能量函數分成權值能量部分項和平滑約束能量部分項。構造圖中的邊分成兩種類型有向邊：列內偏序關系邊和列間限制關系邊。假設每列當中上下相鄰的兩個頂點坐標形式為V（x，y，z），V（x，y，z-1），則列內偏序關系邊可表示為如下形式：

[Edgeintra={V（x，y，z），V（x，y，z-1）|z≥1}] （1）

對于列間限制關系邊，假設在x和y方向的高度差限制分別為Δx和Δy，相鄰列的點構造如下列間限制邊：

[Edgeinter = {V（x，y，z），V（x+1，y，max（z-Δx，0）） V（x，y，z），V（x-1，y，max（z-Δx，0）） V（x，y，z），V（x，y-1，max（z-Δy，0）） V（x，y，z），V（x，y+1，max（z-Δy，0））}] （2）

這兩種關系邊的權值設定為∞，其實際意義在隨后構造s?t圖時該條路徑的流容量無窮大。列間限制為簡單起見采用固定的Δx和Δy，實驗過程令Δx=Δy=3。這兩種限制保證了幀與幀圖像間的限制關系。最終求最優曲面可轉化成求圖的最小閉集。可證明該最小閉集的問題可轉化成求圖的s?t割問題。具體基礎理論內容可參考相應文獻。

1.2 梯度弱化算法

對圖像進行雙邊濾波的醫學圖像去噪[7]，并進行灰度填充預處理去噪之后采用梯度弱化算法。梯度弱化是指將所感興趣之外區域的梯度抑制，降低其他部分對所分割區域的干擾，具體步驟如下。

1.2.1 背景估計

針對每一幀OCT圖像背景和目標閾值分割的閾值通過式（3）確定：

[ti=max0≤k≤255h（Iik）+Δg， i=1，2，…，128] （3）

式中：t表示需要獲取的閾值；i表示體數據的幀數；k表示灰度級；函數[h（k）]計算灰度級k的像素個數；Δg表示容許的灰度差異，可取一個較小的常量，比如Δg=10。通過獲取得到閾值之后對圖像進行閾值分割即可先把背景區域分離出來。

1.2.2 視網膜區域估計

通過視網膜區域估計得到感興趣區域。如上小節所述，得到的目標區域不一定是視網膜區域，也有可能是其他部分的較亮的條紋或者亮斑帶來的干擾。通過取最大連通區域來消除干擾，計算步驟如下：

（1） 對二值圖像連通區域標記。假設二值圖像連通區域個數為n，對應的連通區域分別為[p1，p2，…，pn]，[j∈pipij]統計連通區域p內目標像素個數。二值圖像中兩個目標像素屬于同一個4?連通區域條件是在4?鄰域條件下可達。對每幀二值圖像可通過式（7）計算最大連通分量：

[pmax=maxi=1，2，...，n（j∈pipij）] （4）

（2） 在（1）基礎上，獲取每一列最大連通區域，用來消除視網膜上方玻璃體條帶狀物質帶來的干擾。二值圖像僅保留最大連通分量[pmax]。假設二值圖像第i列連通分量個數為ni，對應的連通區域分別為[ci1 ， ci2，… ，cini]，并且[j∈cimcimj]表示統計連通區域[cim]內目標像素個數。通過式（5）計算每一列的最大連通分量（col表示圖像的列數）：

[cimax=maxm=1，2，...，ni（j∈cimcimj） ， i=1，2，…，col] （5）

在經過以上步驟后得到視網膜目標區域在OCT圖像中的位置。

1.2.3 目標區域腐蝕及梯度弱化

目標區域腐蝕目的是獲得非感興趣區域，并將該區域梯度弱化，主要思想是對得到的視網膜估計的感興趣區域進行形態學腐蝕得到非感興趣區域。假設目標區域的單幀灰度圖像為I，其對應的區域估計后的二值圖像為Ibw。首先通過式（6）對二值圖像進行腐蝕，SE表示腐蝕操作的結構元素，這里采用半徑為3的圓作為結構元素，*為卷積操作。腐蝕操作后的目標區域會收縮變小，從而能夠保證腐蝕后保留出ILM表面附近區域：

[Mdata=Ibw? SE] （6）

[Iout=I*-1-2-1000121?（Mdata）] （7）

將腐蝕之后的目標區域內的梯度通過式（7）去除，其中[Mdata]表示對圖像進行取反操作。去除后保證在輸出圖像中，僅有ILM表面附近的強梯度區域得以保留，最終得到的較大梯度的感興趣區域僅為ILM附近區域。

1.3 ILM表面分割

以弱化后的梯度圖像來計算每個體素的權值，并對整張圖利用3D的圖搜索算法進行操作。采用文獻[9]使用的最大流最小割算法計算最小割。具體步驟如下：

（1） 對圖像4倍下采樣，利用3D圖搜索方法進行分割求解，得到下采樣ILM層分割結果；

（2） 利用對下采樣ILM層分割結果采用文獻[10]使用的薄板樣條函數曲線擬合得到粗分割結果；

（3） 對原始圖像ILM層細分割，在限制范圍內采用3D圖搜索算法，最終得到細分割結果。

2 實驗結果和分析

本文采用了斯坦福大學25對眼睛的SD?OCT圖像進行實驗，圖像體數據分辨率為512×1 024×128。每只眼睛分成OD（右），OS（左）兩組圖像，共50組圖像，利用Matlab腳本編程與C++混編完成算法過程實現。使用ITK?SNAP軟件的snake算法進行半自動分割ILM表面，并利用該軟件本身自帶手動分割功能進行調整后作為標準，與本文算法得出的結果進行對比，同時也對50組圖像直接采用3D圖搜索算法分割ILM表面，給出了本文算法和直接利用3D圖搜索分割ILM層的算法的實驗結果精度對比。

2.1 定性分析

實驗對不同的SD?OCT圖像進行分割，并在圖3中給出了單幀正常圖像的分割結果。在50組數據集中主要包括兩種特殊情況，一是在正常圖像里中央凹區域不明顯，二是分割病變圖像，在ILM層附近主要可能的病變為下積液。圖3給出這兩種特別情況的分割結果。圖3（c） 比圖3（b）分割更加精確，對圖像灰度填充降低了下積液對ILM分割的干擾，不會分割到下積液里去。圖3（f）對中央凹區域梯度不明顯分割效果較精確，原因是結合3D圖搜索本身的列間限制和薄板樣條插值保證分割結果光滑；直接3D圖搜索加薄板樣條插值可能由于中央凹區域較暗，該情況會有可能發生誤分割，圖3（e）給出發生誤分割的情況。

2.2 定量分析

表1給出誤差的平均結果對比（誤差的平均結果不包括誤分割的數據在內），加粗數據為誤分割。實驗采用ITK軟件進行半自動分割得出的ILM表面為基準，分析了本文算法得出ILM表面的位置與基準位置的誤差，并且在圖4中分別給出了OD（右）和OS（左）的對比數據統計結果（由于誤分割均值和標準差較大，為便于觀察對比結果，將誤分割的均值換成了一個相對較小的值8，方差設為0，如圖4所示）。

結合圖4和表1數據，不包含誤分割情況，直接采用3D圖搜索誤差平均結果為OD（右）：（2.68±2.50）μm，OS（左）：（2.94±2.98）μm，本文算法25組OD（右）誤差的平均結果為（2.55±2.35） μm，OS（左）為（2.71±2.73） μm，比直接采用3D圖搜索精度要高，而且避免了誤分割情況的發生。

2.3 分析與討論

（1） 分割精度。由表1可以看出，50組數據中，有47組分割結果比直接利用3D圖搜索算法分割ILM層精度要高，原因是在獲得梯度弱化后的梯度圖像前，算法進一步去除了分割過程中產生的部分干擾，能夠在一定程度上提高分割精度。須注意加粗的異常分割結果，產生原因是直接利用3D圖搜索算法分割ILM層，由于其他層可能梯度較大，而求解過程為獲得全局最優解，在不改變權重的情況下會造成誤分割，如圖5所示。本文算法由于弱化了其他層，避免了誤分割的發生。

（2） 時間復雜度。直接采用3D圖搜索算法分割128幀圖像需要300 s以上的時間；而本文算法所有數據都能在70 s以內分割128幀圖像，平均單張分割不到1 s。原因是由于梯度弱化后在構造s?t圖時由于很多權重為0的頂點不會與源點s和匯點t有邊，從而降低了求解圖割的構圖規模，從而降低了計算的時間復雜度。

（3） 算法魯棒性。能夠精確分割ILM邊界，并對可能出現的中央凹不明顯以及存在下積液孔洞的特殊情況魯棒性較好，如圖5所示。

3 結 語

本文引入梯度弱化的概念分割SDOCT圖像中的ILM層，利用3D圖搜索算法考慮到圖像幀間的互信息的優勢，并通過梯度弱化過程分割ILM表面，降低了其他區域對該表面分割的干擾，提高了分割的精度和圖搜索的速度。但對于圖像體素的權值函數應該結合其他更多的信息。另外，對于3D圖搜索僅采用比較基本的思想，表面限制采用單一的硬限制，今后可考慮將3D圖搜索中的多層同時分割引入梯度弱化思想，考慮構造比較完善的權值函數。

參考文獻

[1] 孫延奎.光學相干層析醫學圖像處理及其應用[J].光學精密工程，2014（4）：1086?1104.

[2] 徐奇，常英，李文彬，等.基于區域生長的OCT圖像分層算法[J].信息技術，2013（4）：136?140.

[3] KAJIC V， ESMAEELPOUR M， POVAZAY B， et al. Automated choroidal segmentation of 1060 nm OCT in healthy and pathologic eyes using a statistical mode [J]. Biomed Opt Express， 2012， 3（1）： 86?103.

[4] 張博聞，田小林，孫延奎.基于改進的數學形態學的OCT圖像快速邊緣提取算法[J].計算機科學，2013（z1）：173?175.

[5] 楊平，彭清，劉維平，等.一種眼底黃斑水腫OCT圖像分割方法[J].生物醫學工程學雜志，2011（5）：1001?1006.

[6] 胡超.基于OCT圖像玻璃疣的自動檢測與分割[D].南京：南京理工大學，2013.

[7] DUFOUR P A， CEKLIC L， ABDILLAHI H. Graph?based multi?surface segmentation of OCT data using trained hard and soft constraints [J]. IEEE Transactions on Medical Imaging， 2013， 32（3）： 531?542.

[8] 張聚，王陳，程蕓.小波與雙邊濾波的醫學超聲圖像去噪[J].中國圖象圖形學報，2014（1）：126?132.

[9] 張鈴.動態網絡上最大流概念及其性質的研究[J].模式識別與人工智能，2013（7）：609?614.

[10] 任偉建，張宏鑫.基于薄板樣條函數的地圖矢量化處理[J].系統仿真技術，2012（3）：187?191.