王培培 孫九愛

1(上海理工大學(xué)醫(yī)療器械與食品學(xué)院，上海 200093) 2(上海健康醫(yī)學(xué)院醫(yī)學(xué)影像學(xué)院，上海 201318)

引言

青光眼是世界上致盲率第二大的眼科疾病，到2010年為止世界范圍內(nèi)大約有6 000萬青光眼患者，到2020年其患病數(shù)量預(yù)計將會繼續(xù)增加2 000萬[1]。青光眼疾病的發(fā)生是由于眼底視神經(jīng)細(xì)胞損傷造成眼底視盤和視杯幾何形狀的改變。視盤也稱視乳頭位于視網(wǎng)膜的遠(yuǎn)端部分，是視神經(jīng)纖維和血管進出眼球的入口，正常情況下由于視神經(jīng)纖維穿過呈現(xiàn)粉紅色。視杯是視盤內(nèi)部無視神經(jīng)纖維通過的部分，呈現(xiàn)蒼白色。青光眼的主要發(fā)病機制為：眼內(nèi)壓升高，篩板受到擠壓導(dǎo)致視神經(jīng)損傷，從而導(dǎo)致視杯變大[2-4]。隨著病程的加長，會出現(xiàn)視覺損傷，甚至是造成不可逆轉(zhuǎn)的失明[5]。在臨床上，一般是通過計算視杯和視盤的垂直距離比來對青光眼進行診斷[6]。因此，視盤、視杯形態(tài)結(jié)構(gòu)的準(zhǔn)確分割和測量對青光眼的診斷是很重要的[7-8]。

目前，在臨床上青光眼的診斷有兩種方式。第一種方式，醫(yī)生依據(jù)眼底照相機拍攝的二維圖像進行判斷，這需要依靠臨床醫(yī)生的行醫(yī)經(jīng)驗，判斷結(jié)果存在一定的主觀性。第二種方式，為了提高診斷精度，同時出現(xiàn)了三維成像設(shè)備，如光學(xué)相干斷層成像(OCT)[9-10]和海德堡視網(wǎng)膜斷層掃描儀(HRT)[11-12]，這可以通過得到眼底的三維圖像對眼底疾病進行診斷。但是，這兩種方式在使用之前都需要對操作人員進行專業(yè)的培訓(xùn)，價格昂貴而且耗時較長[13]。

目前通過眼底照片直接進行三維重建的國內(nèi)外的文獻數(shù)量很少，而且局限性較強。Enrique等主要研究的是眼底局部的三維重建[14]，存在一定的局限性。陳驥等分析了正常眼和屈光異常眼眼底圖像的關(guān)系，建立了從二維平面圖像到三維曲面圖像的映射關(guān)系來對眼底圖像進行三維重建[15]。李超等則是把眼底二維圖像反投影到三維球面上以實現(xiàn)眼底圖像的三維重建[16]。但他們所采用的三維眼底重建方法建立的眼底三維曲面只能描述眼底的總體表面情況，并不能提供更多的三維細(xì)節(jié)信息。

為此，本研究提出了一種眼底圖像三維重建的新方法即從明暗恢復(fù)形狀法(shape from shading，SFS)[17]。其主要思想是，根據(jù)眼底數(shù)字照相機拍攝的二維眼底數(shù)字圖像的灰度值變化恢復(fù)出眼底三維結(jié)構(gòu)形貌。初步實驗結(jié)果表明，該方法可行性較強，能夠提供更多的線索輔助臨床醫(yī)生對青光眼等眼底類疾病進行診斷。

2 方法

2.1 算法實現(xiàn)流程

本算法的實現(xiàn)主要有以下步驟，如圖1所示。第1步，使用數(shù)字眼底照相機拍攝眼底圖像；第2步，對圖像進行預(yù)處理，包括把RGB圖像轉(zhuǎn)化為灰度圖像，進行圖像分割、濾波和白平衡處理；第三步，計算估計眼底照相機拍攝眼底圖像的光源方向；第4步，利用Tsai[18]線性化方法分析求解SFS算法，再結(jié)合第3步求出的光源方向，對眼底圖像進行三維重建；第5步，可視化重建后的視乳頭三維結(jié)構(gòu)形貌圖。

圖1 視乳頭三維重建流程Fig.1 Diagram of reconstructing 3D topographic map of the ONH

2.2 從明暗恢復(fù)形狀法

20世紀(jì)70年代，Horn提出了一種恢復(fù)物體表面三維形貌的新方法——從明暗恢復(fù)形狀法(SFS)[19]。其主要原理是：假設(shè)物體的表面反射模型為朗伯體表面反射模型，計算出圖像表面的光源方向；再根據(jù)二維圖像的明暗(亮度值)變化，經(jīng)過一定的算法，計算出每一點的高度值，恢復(fù)出物體的三維形貌。

一般情況下，當(dāng)光照射到物體表面時，其表面會同時發(fā)生漫反射和鏡面反射。朗伯體表面反射模型是一種理想的情況，其假設(shè)物體表面只有漫反射，沒有鏡面反射，即物體表面的亮度不會隨觀察角度的改變而改變。在使用SFS算法恢復(fù)出物體三維形貌的過程中，需要同時假設(shè)物體表面是光滑的，即其表面函數(shù)具有一階連續(xù)導(dǎo)數(shù)。

SFS是一種病態(tài)問題的逆運算算法，為解決病態(tài)問題，其算法也經(jīng)過了一系列的改進，出現(xiàn)了演化法、最小值法、局部分析法和線性化法?？紤]到實現(xiàn)的效率因素，本研究采用線性化法對SFS算法進行求解，進而實現(xiàn)了眼底視乳頭的三維重建。

2.2.1計算與估計視乳頭光照方向

相比較眼底圖像，視乳頭是一個較小的區(qū)域。視乳頭中粉紅色的部分是血管和視神經(jīng)進入視網(wǎng)膜的一個切入點(視盤)，視盤中沒有視神經(jīng)穿過且呈現(xiàn)灰白色的部分是視杯。在眼底圖像拍攝的過程中，通常以黃斑作為眼底攝影的中心，而視乳頭偏向于鼻側(cè)，因此不能簡單地把視乳頭攝影視為一個同軸眼底光學(xué)系統(tǒng)。所以，在使用SFS方法對視乳頭進行三維重建時，需要預(yù)先估計出光照方向[20]。

圖2 反射模型Fig.2 Surface reflection model

圖2所示為視乳頭的表面反射模型。以二維圖像上的某一點為例，建立直角坐標(biāo)系，其中z軸表示物體的高度方向，x軸的值表示圖像上某一點像素點的行數(shù)，y軸的值代表圖像上某一點像素點的列數(shù)。假設(shè)α角為x軸與表面法向量在xy平面的投影的夾角，β表示z軸與表面法向量的夾角，τ表示光源方向在xy平面的投影與x軸的夾角，σ表示光源方向與z軸的夾角。其中，α,τ∈(0,2π)，β,σ∈(0,π/2)。

根據(jù)笛卡爾坐標(biāo)系和球面坐標(biāo)系的關(guān)系，視乳頭的表面法向量，光源方向和輻射方程可以分別表示為

n=[sinβcosα,sinβsinα,cosβ]T

(1)

I=[sinσcosτ,sinσsinτ,cosσ]T

(2)

E(α,β)=(sinβcosαsinσcosτ,

sinβcosαsinσsinτ,cosβcosσ)

(3)

β與cosβ呈線性變化，α服從均勻分布。σ和τ可以通過對亮度方程和亮度的平方方程求解積分得到，即

f(α,β)=cosβ/2π

(4)

式中，μ1和μ2分別表示所給圖像的亮度值和亮度值的平方。

利用式(5)、(6)可求出σ和τ的表達式，有

(7)

(8)

利用式(7)、(8)即可求出其角度值，把其角度值代入式(2)中，即可得到光源方向。

2.2.2SFS算法實現(xiàn)

假設(shè)視乳頭區(qū)域表面的反射模型為朗伯體反射模型，則圖像的亮度方程為

(9)

或

(10)

其中，令I(lǐng)=(n01,n02,n03)和I=(-p0,-q0,1)表示光源方向，兩者的關(guān)系為

(11)

(p,q)為圖像上任意一點的梯度值，任意一點的高度值為z=z(x，y)，表面法向量為n=(n1,n2,n3)。高度值和梯度值之間的關(guān)系為p=?z/?x，q=?z/?y。表面法向量和梯度值是使用SFS算法進行三維重建的兩個重要的參數(shù)值。假設(shè)一幅圖像有M×N個像素值，根據(jù)SFS算法其有M×N個方程，但是卻有2×M×N個待求量。為解決這個問題，本研究采用Tsai的離散線性化方法對方程進行求解。

根據(jù)Tsai的離散化原理，令z(i,j)表示圖像上第(i,j)點的高度值，則圖像的表面梯度可以表示為

(12)

假設(shè)一幅圖像的行數(shù)和列數(shù)分別為M、N，則i和j的取值范圍為:i=0, 1, 2，…，M-1;j=0, 1, 2,…,N-1。則亮度方程的差分方程為

(13)

再利用泰勒公式改寫式(13)，得到

(14)

(15)

其中

(16)

3 實驗結(jié)果

為得到局部視杯圖像，本研究對RIM-ONE眼底圖像數(shù)據(jù)庫的圖像[21](包含正常眼和疑似青光眼圖像)進行預(yù)處理，利用自動閾值分割和區(qū)域生長法對眼底圖像進行分割，得到了如圖3所示的包含視杯的局部眼底圖像。

圖3 眼底RGB圖像。(a)和(b)是正常眼底視杯圖像；(c)和(d)是疑似青光眼患者眼底圖像Fig.3 The RGB fundus image.(a)and (b) are normal fundus images; (c) and (d) are suspicious glaucoma fundus images

利用Matlab軟件對圖像進行白平衡處理，消除光線對圖像的影響，然后把RGB圖像轉(zhuǎn)化為灰度圖像，接著利用圖像的亮度(灰度值信息)求出每幅視乳頭圖像的光源方向。通過對式(5)、(6)進行求解，計算出4幅圖像的光源方向為

Ia=(0.543 1,-0.023 9,0.839 3)

Ib=(0.339 6,0.523 8,0.839 3)

Ic=(0.566 7,-0.266 1,0.792 3)

Id=(0.616 9,-0.030 3,0.786 4)。

然后，利用上述所求的光源方向結(jié)合SFS算法，對圖像的每一個像素點的灰度值進行處理，得到圖像三維數(shù)據(jù)，對眼底表面形貌進行三維重建，即可得到如圖4所示的重建后的視乳頭三維形貌圖。

圖4 重建后的三維眼底圖。(a)和(b)是正常眼底視杯圖像；(c)和(d)是疑似青光眼患者眼底圖像Fig.4 Reconstructed three dimensional fundus image.(a)and (b) The normal fundus images; (c) and (d) The suspicious glaucoma fundus image

在圖4中，(a)和(b)是正常眼底圖像的視杯圖形，(c)和(d)是青光眼患者的眼底圖像。從圖4中可以很清晰地分辨出(c)和(d)兩幅青光眼眼底形貌圖中視盤和視杯的邊界信息。計算每幅眼底圖像中視盤和視杯的垂直直徑比可得出，圖4中的眼底圖像(a)～(d)的杯盤比比值分別為0.46、0.35、0.611、0.614。由此亦可判斷出，圖4中的前兩幅圖像為正常眼的視乳頭三維重建結(jié)果，后兩幅為青光眼患者的視乳頭三維形貌圖。此結(jié)果與事實相符合。

因此，視乳頭形貌的三維重建，可以更加方便、準(zhǔn)確地計算視杯和視盤的比值，對青光眼地診斷和輔助治療提供幫助。

4 討論和結(jié)論

本研究主要提出了一種獲取眼底結(jié)構(gòu)三維形貌的方法，為眼底疾病的診斷和跟蹤治療提供一種簡單、經(jīng)濟方便的手段。而且，從圖4的三維重建結(jié)果可以看出，利用SFS算法對眼底局部圖像進行三維重建的效果較好，能夠很清晰地提供視盤和視杯的邊界信息，更加準(zhǔn)確地計算出視杯和視盤比值，為青光眼的診斷和跟蹤治療提供更加有效便利的方式。

在本研究中，利用閾值分割和區(qū)域生長法，對數(shù)字眼底照相機獲得的眼底圖像進行分割，得到視乳頭區(qū)域圖像，如圖3所示。再利用白平衡消除光線影響，然后經(jīng)圖像濾波和灰度變換得到眼底灰度圖像。接著利用圖像亮度方程和亮度值得到視乳頭圖像光源方向。但是，使用此方法雖能得到眼底圖像的光源方向，但其結(jié)果受到圖像預(yù)處理的影響。最后利用SFS算法計算得出視乳頭圖像的三維高度值如圖4中的圖像所示。

臨床上，當(dāng)視杯視盤的比值大于等于0.6時[6]，即可診斷為青光眼。從本研究的實驗結(jié)果中得出的視乳頭中視杯和視盤的比值可以發(fā)現(xiàn)，利用本研究的算法得出的青光眼患者的視乳頭三維重建形貌圖形中的杯盤比大于0.6，正常眼的杯盤比小于0.6，與臨床上的診斷結(jié)果相吻合。因此，本算法系統(tǒng)可以輔助青光眼的診斷。

但是，目前的研究只對與青光眼有關(guān)的視盤部分進行了三維重建。下一步將擴展算法，以便能夠?qū)ν暾鄣讏D像進行重建，為更多眼底疾病的診斷和跟蹤治療提供可能。

[1] Quigley HA, Broman AT. The number of people with glaucoma worldwide in 2010 and 2020 [J]. Digest of the World Core Medical Journals, 2006,90(3):262-267.

[2] 黨鴻, 辛?xí)匀? 青光眼視神經(jīng)損傷機制的研究進展 [J]. 眼科新進展, 2016, 36(7):680-683.

[3] Septiarini A, Harjoko A. Automatic glaucoma detection based on the type of features used: A review [J]. Journal of Theoretical & Applied Information Technology, 2015,2872(3): 366-375.

[4] Kilintzis V, Chouvarda I, Topouzis F, et al. Symbolic representation of optic nerve head cup presents difference in patients with Primary Open Angle Glaucoma[C]// IEEE International Conference on Information Technology and Applications in Biomedicine. Corfu:IEEE, 2011:1-4.

[5] Heijl A, Aspberg J, Bengtsson B. The effect of different criteria on the number of patients blind from open-angle glaucoma[J]. BMC Ophthalmology, 2011, 11(1):1-6.

[6] Almazroa A, Burman R, Raahemifar K, et al. Optic disc and optic cup segmentation methodologies for glaucoma image detection: A survey [J]. Journal of Ophthalmology, 2015, 2015(7):1-28

[7] Mary MCVS, Rajsingh EB, Jacob JKK, et al. An empirical study on optic disc segmentation using an active contour model [J]. Biomedical Signal Processing & Control, 2015, 18:19-29.

[8] Nawaldgi S. Review of automated glaucoma detection techniques [C]//International Conference on Wireless Communications. Signal Processing and Networking. Boston:IEEE, 2016:1435-1438.

[9] Huang D, Swanson EA, Lin CP et al. Optical coherence tomography[J]. Science, 1991,254(5035):1178-1181.

[10] 郭慧敏. OCT檢測視網(wǎng)膜神經(jīng)纖維層厚度在青光眼診斷中的應(yīng)用進展 [J]. 醫(yī)學(xué)綜述, 2013, 19(7):1281-1283.

[11] Kamal DS,Garwayheath DF,Hitchings RA,et al. Use of sequential Heidelberg retina tomograph images to identify changes at the optic disc in ocular hypertensive patients at risk of developing glaucoma[J]. British Journal of Ophthalmology, 2000,84 (9): 993-998.

[12] 邢業(yè)嬌, 王大博, 紀(jì)珍,等. 海德堡OCT測量后極部視網(wǎng)膜厚度對青光眼診斷價值 [J]. 青島大學(xué)醫(yī)學(xué)院學(xué)報.2013, 2013(1):38-40.

[13] Yin F, Wong DWK, Quan Y, Ai PY, Tan NM, Gopalakrishnan K, et al. A cloud-based system for automatic glaucoma screening [C]//Annual International Coference of the IEEE. Milan:IEEE, 2015: 1596-1599.

[14] Corona E, Mitra S, Wilson M, et al. Digital stereo image analyzer for generating automated 3-D measures of optic disc deformation in glaucoma[J]. IEEE Transactions on Medical Imaging, 2002, 21 (10): 1244-1253.

[15] 陳驥, 彭承琳. 眼底圖像的三維重建[J]. 生物醫(yī)學(xué)工程雜志, 2008,25(1):177-181.

[16] 李超, 梁斌, 陳武凡, 等. 由二維眼底正投影圖像向三維曲面逆投影成像的重建算法[J]. 中國生物醫(yī)學(xué)工程學(xué)報, 2002,21(4):346-350.

[17] Horn BKP. Height and gradient from shading[J]. International Journal of Computer Vision, 1990,5(1):37-75.

[18] Tsai PS, Shah M, editors. A fast linear shape from shading [C]// Proceedings 1992 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Champaign:IEEE, 1992:734-736.

[19] Horn BKP. Shape from Shading: A Method for Obtaining the Shape of a Smooth Opaque Object from One View[M]. Cambridge: MIT Press,1970.

[20] Elhabian SY. Hands on shape from shading [J]. Technical Report, 2008.

[21] Fumero F, Alayon S, Sanchez JL, et al. RIM-ONE: An open retinal image database for optic nerve evaluation[C]// Proceedings of the 2011 24th International Symposium on Computer-Based Medical Systems. Washington DC: IEEE, 2011:1-6.