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

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

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

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

0 引 言

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

近年來針對視網(wǎng)膜分層的方法主要包含如下幾類。文獻[2]提出基于區(qū)域增長的分割方法，該方法不能保證所求表面全局最優(yōu)；Kajic等人給出了一種基于圖像紋理和形狀統(tǒng)計信息的邊緣特征提取分割算法，醫(yī)學圖像的紋理特征提取相對比較困難[3]； 文獻[4]采用結(jié)合形態(tài)學的邊緣提取方法，該方法針對部分邊緣不明顯的醫(yī)學圖像效果不甚好，沒有利用幀間信息；文獻[5]給出了基于水平集方法分割OCT圖像，其對應的符號距離函數(shù)的計算量較大。以上幾類分割方法均以O(shè)CT圖像的一個B?Scan（一幀切片圖像）為單位進行圖像分割，沒有充分考慮連續(xù)圖像幀與幀之間相對連續(xù)關(guān)系；目前基于圖論的分割方法可分成2D和3D兩類。2D的圖論算法的圖像分割單位針對單幀圖像。文獻[6]采用圖搜索的方法對OCT視網(wǎng)膜圖像進行層分割，利用Dijkstra算法搜索出最短路徑，逐步區(qū)域限制分割出各個視網(wǎng)膜層，其針對一些邊緣不明顯圖像也需要做部分分割后處理操作?；?D圖論的方法在近年來得到廣泛的應用。3D圖搜索采用構(gòu)圖的方法，給出了多表面同時分割的理念并引入3D子圖間限制，最終轉(zhuǎn)換成求解一個4D結(jié)構(gòu)圖的最優(yōu)解問題，目前Dufour等人在3D圖搜索基礎(chǔ)上又給出了先驗信息模型的概念，并結(jié)合先驗信息模型給出了軟限制構(gòu)圖的思想，但依然需要層與層之間約束構(gòu)圖同時分割求解，不能做到簡單，單獨分割某層[7]。如圖1所示。其中表面標記從上到下分別為內(nèi)界膜（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圖搜索算法求解較慢，并且難以單獨分割某一層。本文在結(jié)合3D圖搜索方法基礎(chǔ)上，引入一種梯度弱化的方法，并將該方法應用到ILM表面的分割問題上。該方法優(yōu)勢是，第一，能夠單獨分割I(lǐng)LM層。直接單獨分割I(lǐng)LM層可能會產(chǎn)生誤分割，常采用3層同時分割才可分出ILM層，當僅需分割I(lǐng)LM層時不甚合理；第二，利用梯度弱化后結(jié)合3D圖搜索分割I(lǐng)LM層會加快圖搜索的速度，原因是僅保持ILM表面附近較大梯度，其他區(qū)域梯度較小且相對一致，頂點權(quán)重w的計算方式導致構(gòu)圖過程很多頂點權(quán)重為0，在將3D圖搜索轉(zhuǎn)化成求解最小割問題后，這些點會排除在s?t圖以外，大大減少s?t圖中的邊數(shù)，加快求解最小割速度。

1 方 法

本節(jié)給出梯度弱化方法的概念及算法實現(xiàn)。算法實現(xiàn)流程如圖2所示。

1.1 3D圖搜索算法模型

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

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

對于列間限制關(guān)系邊，假設(shè)在x和y方向的高度差限制分別為Δx和Δy，相鄰列的點構(gòu)造如下列間限制邊：

[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）

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

1.2 梯度弱化算法

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

1.2.1 背景估計

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

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

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

1.2.2 視網(wǎng)膜區(qū)域估計

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

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

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

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

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

在經(jīng)過以上步驟后得到視網(wǎng)膜目標區(qū)域在OCT圖像中的位置。

1.2.3 目標區(qū)域腐蝕及梯度弱化

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

[Mdata=Ibw? SE] （6）

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

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

1.3 ILM表面分割

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

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

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

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

2 實驗結(jié)果和分析

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

2.1 定性分析

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

2.2 定量分析

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

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

2.3 分析與討論

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

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

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

3 結(jié) 語

本文引入梯度弱化的概念分割SDOCT圖像中的ILM層，利用3D圖搜索算法考慮到圖像幀間的互信息的優(yōu)勢，并通過梯度弱化過程分割I(lǐng)LM表面，降低了其他區(qū)域?qū)υ摫砻娣指畹母蓴_，提高了分割的精度和圖搜索的速度。但對于圖像體素的權(quán)值函數(shù)應該結(jié)合其他更多的信息。另外，對于3D圖搜索僅采用比較基本的思想，表面限制采用單一的硬限制，今后可考慮將3D圖搜索中的多層同時分割引入梯度弱化思想，考慮構(gòu)造比較完善的權(quán)值函數(shù)。

參考文獻

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

[2] 徐奇，常英，李文彬，等.基于區(qū)域生長的OCT圖像分層算法[J].信息技術(shù)，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] 張博聞，田小林，孫延奎.基于改進的數(shù)學形態(tài)學的OCT圖像快速邊緣提取算法[J].計算機科學，2013（z1）：173?175.

[5] 楊平，彭清，劉維平，等.一種眼底黃斑水腫OCT圖像分割方法[J].生物醫(yī)學工程學雜志，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] 張聚，王陳，程蕓.小波與雙邊濾波的醫(yī)學超聲圖像去噪[J].中國圖象圖形學報，2014（1）：126?132.

[9] 張鈴.動態(tài)網(wǎng)絡(luò)上最大流概念及其性質(zhì)的研究[J].模式識別與人工智能，2013（7）：609?614.

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