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

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

中圖分類號(hào)： TN919?34； TP391 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2015）04?0106?05

0 引 言

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

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

1 方 法

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

1.1 3D圖搜索算法模型

3D圖搜索算法能量函數(shù)分成權(quán)值能量部分項(xiàng)和平滑約束能量部分項(xiàng)。構(gòu)造圖中的邊分成兩種類型有向邊：列內(nèi)偏序關(guān)系邊和列間限制關(guān)系邊。假設(shè)每列當(dāng)中上下相鄰的兩個(gè)頂點(diǎn)坐標(biāo)形式為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）

對(duì)于列間限制關(guān)系邊，假設(shè)在x和y方向的高度差限制分別為Δx和Δy，相鄰列的點(diǎn)構(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è)定為∞，其實(shí)際意義在隨后構(gòu)造s?t圖時(shí)該條路徑的流容量無窮大。列間限制為簡單起見采用固定的Δx和Δy，實(shí)驗(yàn)過程令Δx=Δy=3。這兩種限制保證了幀與幀圖像間的限制關(guān)系。最終求最優(yōu)曲面可轉(zhuǎn)化成求圖的最小閉集。可證明該最小閉集的問題可轉(zhuǎn)化成求圖的s?t割問題。具體基礎(chǔ)理論內(nèi)容可參考相應(yīng)文獻(xiàn)。

1.2 梯度弱化算法

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

1.2.1 背景估計(jì)

針對(duì)每一幀OCT圖像背景和目標(biāo)閾值分割的閾值通過式（3）確定：

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

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

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

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

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

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

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

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

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

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

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

[Mdata=Ibw? SE] （6）

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

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

1.3 ILM表面分割

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

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

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

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

2 實(shí)驗(yàn)結(jié)果和分析

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

2.1 定性分析

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

2.2 定量分析

表1給出誤差的平均結(jié)果對(duì)比（誤差的平均結(jié)果不包括誤分割的數(shù)據(jù)在內(nèi)），加粗?jǐn)?shù)據(jù)為誤分割。實(shí)驗(yàn)采用ITK軟件進(jìn)行半自動(dòng)分割得出的ILM表面為基準(zhǔn)，分析了本文算法得出ILM表面的位置與基準(zhǔn)位置的誤差，并且在圖4中分別給出了OD（右）和OS（左）的對(duì)比數(shù)據(jù)統(tǒng)計(jì)結(jié)果（由于誤分割均值和標(biāo)準(zhǔn)差較大，為便于觀察對(duì)比結(jié)果，將誤分割的均值換成了一個(gè)相對(duì)較小的值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層精度要高，原因是在獲得梯度弱化后的梯度圖像前，算法進(jìn)一步去除了分割過程中產(chǎn)生的部分干擾，能夠在一定程度上提高分割精度。須注意加粗的異常分割結(jié)果，產(chǎn)生原因是直接利用3D圖搜索算法分割I(lǐng)LM層，由于其他層可能梯度較大，而求解過程為獲得全局最優(yōu)解，在不改變權(quán)重的情況下會(huì)造成誤分割，如圖5所示。本文算法由于弱化了其他層，避免了誤分割的發(fā)生。

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

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

3 結(jié) 語

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

參考文獻(xiàn)

[1] 孫延奎.光學(xué)相干層析醫(yī)學(xué)圖像處理及其應(yīng)用[J].光學(xué)精密工程，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] 張博聞，田小林，孫延奎.基于改進(jìn)的數(shù)學(xué)形態(tài)學(xué)的OCT圖像快速邊緣提取算法[J].計(jì)算機(jī)科學(xué)，2013（z1）：173?175.

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

[6] 胡超.基于OCT圖像玻璃疣的自動(dòng)檢測與分割[D].南京：南京理工大學(xué)，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ī)學(xué)超聲圖像去噪[J].中國圖象圖形學(xué)報(bào)，2014（1）：126?132.

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

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