基于圖像處理技術的瞳孔和角膜反射中心提取算法

王錦榕，袁學海，劉增良

(1.大連理工大學控制科學與工程學院，遼寧大連116024;2.中國人民解放軍國防大學信息作戰(zhàn)與指揮訓練教研部，北京100091)

人總是通過五官從周圍世界獲取信息，并根據(jù)這些信息指導自己的行為，而其中80%～90%的信息是由視覺來獲得的，因此，對如何利用眼睛獲取視覺信息進行研究是至關重要的.眼動測量一直是國內外有關學者關注與研究的熱點課題.近年來，隨著傳感技術和微電子技術的飛速發(fā)展，眼動測量技術得到了很大的發(fā)展，測量眼球運動的方法已日趨成熟，在工業(yè)控制、機器人學和臨床醫(yī)學上有著廣泛的應用前景.目前主要的測量方法有眼電圖法[1]、角膜反射法[2]、電磁線圈法[3]、紅外光電反射法[4]、紅外 TV 法[5]、角膜反射法[6]、瞳孔－角膜跟蹤法[7]、Purkinje影像追蹤法[7]等.在以上各種方法中，瞳孔－角膜跟蹤法測量準確、誤差小，且對人無干擾，是目前較為理想的一種眼動測量方法.

文獻[8]中采用星射線方法獲取瞳孔邊界點，但在眼瞼遮擋部分瞳孔的情況下，易受到眼瞼和睫毛的干擾，很難獲取準確的瞳孔邊界點.基于最小二乘法和隨機原理的橢圓擬合方法[9-10]，通過對48個瞳孔邊界特征點進行橢圓擬合，考慮到算法的準確性和實時性，執(zhí)行的次數(shù)一般取100～150之間，該算法的執(zhí)行時間約為8 ms.為了進一步減小瞳孔橢圓擬合時間，以及在眼瞼遮擋部分瞳孔的情況下實現(xiàn)瞳孔中心的準確定位，本文提出了一種基于圖像處理技術的瞳孔和角膜反射中心提取算法.

1 角膜反射中心和瞳孔中心的提取

在紅外光源下，用紅外濾光片濾去可見光，通過攝像機實時采集人眼圖像視頻中的人眼圖像，見圖1(a).圖像中有一個很亮的反射光斑，即角膜反射光斑.瞳孔的形狀近似為橢圓，顏色較暗，其外圍虹膜顏色較亮，瞳孔與虹膜的邊界清晰可見.角膜反射中心和瞳孔中心的提取過程可分為6個步驟:圖像預處理、獲取圖像處理區(qū)域、提取角膜反射中心、求取自適應最佳閾值確定瞳孔位置和大小、用梯度法獲取瞳孔輪廓特征點、使用橢圓擬合瞳孔輪廓確定瞳孔中心.如圖1所示.

圖1 圖像處理結果Fig.1 The results of image processing

1.1 圖像預處理

圖像在拍攝和傳送過程中會不可避免地受到噪聲的干擾，因此先要對圖像進行平滑處理以抑制噪聲，改善圖像質量.常用的平滑處理方法有鄰域平均法、中值濾波、高斯濾波等.由于高斯濾波能有效地去除服從正態(tài)分布的噪聲，且不會被無用的高頻信號污染，同時保留了大部分所需要的信號，故采用高斯濾波對人眼圖像進行平滑處理.

1.2 獲取圖像處理區(qū)域

在圖像處理過程中，主要任務是提取角膜反射中心和瞳孔中心，因此只需要對圖像中眼睛所在的區(qū)域進行處理即可，從而提高處理速度.由于人眼圖像中，瞳孔顏色較暗，即灰度值較低，因此可以采用具有自適應于每幀圖像變化的閾值分割法，來獲取瞳孔所在的區(qū)域.

一般把包含某個目標物體的矩形外框稱為該物體所在的矩形區(qū)域，它由矩形的左上角坐標和右下角坐標表示.獲取圖像處理區(qū)域的算法如下.

1)取一個較低的閾值t，求出灰度值低于t的面積s和所在的矩形區(qū)域R.

2)當s小于某個設定門限閾值時，t加1，再求出 s和 R.

3)重復執(zhí)行步驟2)，當t達到某個設定值，且s小于設定門限值時，則認為眼睛閉著，退出循環(huán);反之，則繼續(xù)循環(huán)，直到s達到設定門限時求出R，并退出.

將R各邊擴大80個像數(shù)，使瞳孔和角膜反射光斑在R內，令R為圖像處理區(qū)域，處理結果見圖1(b).

1.3 提取角膜反射中心

由圖1(a)可知，角膜反射光斑最大、中心最亮，從中心由里往外亮度逐漸變暗，其中還有一些較小的雜質光斑.當取較大的圖像二值化閾值時，角膜反射光斑的面積較小，幾乎沒有雜質點，但提取的角膜反射中心誤差較大.隨著二值化閾值的減小，角膜反射光斑的面積逐漸變大，同時引入的雜質點也逐漸增加.因此，要使提取的角膜反射中心誤差減小，必須使角膜反射光斑面積足夠大，且雜質點很少.

由于角膜反射光斑中心最亮，因此可以用高低2次閾值提取角膜反射中心，具體算法如下.

1)首先取一個較大的圖像二值化閾值，求取粗略的角膜反射中心.

2)以這個角膜反射中心為圓心，r為半徑的圓形區(qū)域作為圖像處理區(qū)域，其中r的取值要使得角膜反射光斑在圓形區(qū)域內.

3)在圖像處理區(qū)域內，取一個較小的二值化閾值，求取角膜反射中心坐標(xc，yc).

此方法增大了角膜反射光斑面積的同時，消除了雜質光斑的干擾，既提高了角膜反射中心的提取精度，同時又減少了處理時間，處理結果見圖1(c).

1.4 求取自適應閾值確定瞳孔位置和大小

由于瞳孔較暗，灰度值比較小，因此可以采用具有自適應于每幀圖像變化的閾值來提取瞳孔.當二值化閾值很小時，只有灰度值小于閾值的點可以通過，則瞳孔面積和瞳孔所在的矩形區(qū)域很小，是整個瞳孔的一小部分.隨著閾值的增加，瞳孔面積和瞳孔所在的矩形區(qū)域面積逐漸加大.由于瞳孔分布均勻，灰度值大小差不多，故瞳孔面積增大的速度比矩形區(qū)域快，即瞳孔面積與矩形區(qū)域面積的比值逐漸增大.當閾值達到一定值時，即擴大到瞳孔的邊緣時，瞳孔面積增加速度逐漸變慢，而矩形區(qū)域面積則會繼續(xù)增加，若繼續(xù)增加閾值，則瞳孔面積與矩形區(qū)域面積的比值逐漸減小.該比值可以表示為

式中:S0表示瞳孔面積，S表示矩形區(qū)域面積.隨著二值化閾值的逐漸增大，p是一個先增后減的比率函數(shù)，此函數(shù)極點所對應的閾值就是提取瞳孔的最佳二值化閾值.通過這個最佳閾值就可以求出瞳孔的位置和大小.

由于人的眼瞼張合度的不同，會出現(xiàn)眼睫毛遮住瞳孔的現(xiàn)象，造成瞳孔中心與矩形區(qū)域中心位置相差很大.在處理過程中，取一個大于瞳孔實際大小的設定值，當瞳孔所在的矩形區(qū)域的長度或寬度大于該設定值時，則認為眼睫毛遮住了瞳孔.由于瞳孔中心與矩形區(qū)域中心位置相差很大，因此可以通過瞳孔中心到矩形區(qū)域各邊的距離來重新設置圖像處理區(qū)域，從而消除眼睫毛的影響，見圖2.

在圖像處理區(qū)域中，求取自適應最佳閾值的具體步驟如下.

1)設置初始二值化閾值t=25，求出初始比率p和矩形區(qū)域R，其中矩形區(qū)域由矩形的左上角坐標(x10，y10)和右下角坐標(x20，y20)表示.

2)當t＜55時，t加1，求出新的比率 p1和矩形區(qū)域R1，R1由矩形的左上角坐標(x11，y11)和右下角坐標(x22，y22)表示，否則退出.

3)當 p－p1＞0.05且 p＞0.6，則此時的 t為最佳閾值，退出;反之，則繼續(xù)執(zhí)行.

4)求R1的長度a和寬度b，若滿足a＞110或b＞100，則繼續(xù)執(zhí)行;反之，令 p=p1，R=R1返回到步驟2)執(zhí)行.

5)令 R0=R1，R0由矩形的左上角坐標(x1，y1)和右下角坐標(x2，y2)表示.求出瞳孔中心坐標(x，y)，再求瞳孔中心到矩形區(qū)域各邊的距離，即lx=x － x1，rx=x2－ x，ty=y －y1，by=y2－ y.

6)若滿足 a＞110，則當 lx＞rx時，令 x1=x－rx，x2=x2+10，y1=y1－10，y2=y2+10;當 lx＜rx時，令x2=x+lx，x1=x1－10，y1=y1－10，y2=y2+10.

7)若滿足b＞100，則當 ty＜by時，令 y1=y1－10，y2=y+ty;當 ty＞ by時，令 y1=y－ by，y2=y2+10.

8)令R0為圖像處理區(qū)域，t減2，求出p和 R，返回到步驟2)執(zhí)行.

通過求出的自適應最佳閾值確定瞳孔位置和大小，該方法具有很強的魯棒性，其處理結果見圖1(d).

圖2 設置處理區(qū)域Fig.2 Set the processing area

1.5 用梯度法提取瞳孔輪廓特征點

在人眼圖像中，瞳孔和虹膜具有明顯的顏色差異，瞳孔的灰度值較低，虹膜的灰度值較高.因此可以通過灰度的變化率和變化方向來衡量瞳孔與虹膜邊界的灰度變化，這里分別用梯度的幅值和方向表示.對于某一方向上一定范圍內的連續(xù)像素點集合，可以用特定的算子求取任意像素點在該方向上的梯度值，梯度值最大的像素點為邊界點[11].

通過自適應最佳閾值法確定的瞳孔中心坐標為(xp，yp)，瞳孔在x方向和y方向上的半徑分別為rx和ry，這些值具有很高的精度和穩(wěn)定性.令

因此，瞳孔輪廓特征點一般在以(xp，yp)為圓心，半徑為r的圓附近.用梯度法提取瞳孔輪廓特征點的具體步驟如下.

1)以(xp，yp)為圓心，半徑為r0的圓環(huán)區(qū)域作為圖像處理區(qū)域，見圖3，其中r－7≤r0≤r+7.

圖3 圓環(huán)區(qū)域Fig.3 The ring area

2)如圖3，在處理區(qū)域內，每一組角度θ和半徑ri都可以在圓上確定一個相應的像素點pi.其中ri=r－7+i(i=0，1，…，14)，θ為眼圖坐標系中，以(xp，yp)為出發(fā)點的射線與橫軸之間的夾角，0≤θ＜2π.這樣可以求出像素點pi在圖像中的坐標為(xi，yi)，其中，

3)將像數(shù)點 pi(i=1，2，…，14)與圓心相連，在連線上取pi的相鄰像數(shù)點pi－1，求取pi在θ方向上的梯度值Ti.設像數(shù)點pi的灰度值為Li，則

令梯度值最大的像數(shù)點為θ方向上瞳孔與虹膜的邊界點.

4)取一組不同的θ，執(zhí)行步驟2)和3)就可以求出一組瞳孔輪廓的特征點.

通過此方法求出的特征點是一組準確的邊界點，其處理結果見圖1(e).

1.6 橢圓擬合瞳孔輪廓

提取的瞳孔輪廓特征點是一組離散數(shù)據(jù)，可以通過橢圓擬合來確定瞳孔中心.基于最小二乘法和隨機原理擬合橢圓的方法具有很高的精度，但需要一定的執(zhí)行時間且魯棒性較差，并且在實時系統(tǒng)中，容易出現(xiàn)干擾數(shù)據(jù).由于用梯度法求取的特征點是一組精確的邊界點，故可以采用基于幾何參數(shù)擬合橢圓的方法，橢圓的標準方程為

式中:橢圓中心坐標為(x0，y0)，長半軸為a，短半軸為b，長軸的轉角為θ.在人眼圖像中，令θ=0，具體的橢圓擬合步驟如下.

1)求取瞳孔輪廓特征點的幾何中心(xp，yp)，以(xp，yp)為橢圓中心，a=rx，b=ry擬合橢圓.

2)統(tǒng)計所有特征點與橢圓的距離小于設定值ε的特征點數(shù)目n，當n大于設定值時，則擬合成功，執(zhí)行步驟5);否則繼續(xù)執(zhí)行步驟3).

3)以(xp，yp)為橢圓中心，nmax=0，令 a 的初始值為r－7，終止值為 r+7，步長為1，b的初始值為r－7，終止值為r+7，步長為1，進行雙重循環(huán).求取每一組a和b對應的n值，當nmax＜n時，令nmax=n，并保存這一組的a和b.

4)當nmax大于設定值時，則擬合成功，保存的a和b就是最佳的長半軸和短半軸，執(zhí)行步驟5);反之，擬合失敗，并退出.

5)除去與橢圓的距離大于設定值ε的特征點，重新求取特征點的幾何中心為瞳孔中心.

對150張圖像進行瞳孔橢圓擬合，基于最小二乘法和隨機原理擬合橢圓方法和基于幾何參數(shù)擬合橢圓方法的比較結果如表1.從表1中可以看出，該方法具有很高的擬合精度，且執(zhí)行時間短，其處理結果見圖1(f).

表1 2種擬合算法的結果比較Table 1 Comparison of the results of two algorithms

2 實驗結果與分析

瞳孔和角膜反射中心的提取算法在以VC++6.0為平臺的P4 2.4GHz主頻、512MB內存的PC機上實現(xiàn).用該算法對人眼圖像進行處理，得到準確的瞳孔中心和角膜反射中心以及最優(yōu)的瞳孔橢圓擬合，見圖1(f).在眼瞼遮擋部分瞳孔的情況下，該算法也能擬合出正確的瞳孔輪廓，見圖1(g)和(h).此算法的平均運行時間為15 ms，完全滿足實時系統(tǒng)要求.因此可以通過此算法對從攝像機實時采集的人眼圖像進行處理，獲取瞳孔中心坐標(xp，yp)和角膜反射中心坐標(xc，yc).瞳孔中心相對于角膜反射中心的偏移量稱為角膜－瞳孔坐標(xcp，ycp)，即

眼睛注視運動和平滑跟蹤運動[7]所對應的角膜－瞳孔位置坐標和實時曲線見圖4.圖4(a)和(b)為眼睛注視前方某一目標物體所得到的角膜－瞳孔位置坐標和實時曲線;圖4(c)和(d)為眼睛圍繞著液晶顯示屏外框做平滑跟蹤運動所得到的角膜－瞳孔位置坐標和實時曲線.實驗表明，在眼瞼遮住部分瞳孔的情況下新算法實現(xiàn)了瞳孔中心的準確提取，且滿足實時性要求.由角膜－瞳孔坐標的實時曲線可知，提取的瞳孔中心和角膜反射中心在實時系統(tǒng)中具有較強的抗干擾能力.

圖4 注視運動和平滑跟蹤運動Fig.4 The fixation and smooth pursuit

3 結束語

瞳孔中心是確定視線方向最主要的因素，能否提取穩(wěn)定、準確的瞳孔中心直接影響視線方向的準確性.本算法采用自適應閾值法確定瞳孔的位置，大大減少了處理的數(shù)據(jù)量，然后用梯度法提取瞳孔輪廓特征點，避免了因眼瞼、眼睫毛干擾產(chǎn)生的虛假特征點，最后基于幾何參數(shù)法的橢圓擬合提取了準確的瞳孔中心，并有效地減少了橢圓擬合時間.后續(xù)將通過角膜－瞳孔坐標確定視線方向，進一步研究眼控鼠標.

[1]ROBINSON D A.A method of measuring eye movements using a scleral search coil in a magnetic field[J].IEEE Transactions on Biomedical Engineering，1963，BME-10:137-145.

[2]SHUPERT C，F(xiàn)UCHS A F.Development of conjugate human eye movements[J].Vision Research，1988，28(5):585-596.

[3]GAUTHIER G M.Vertical and horizontal eye movement recording in the unrestrained cat[J].IEEE Transactions on Biomedical Engineering，1975，BME-22(3):202-207.

[4]JONES R.Two dimensional eye movement recording using a photoelectric matrix method[J].Vision Research，1973，13(2):425-431.

[5]王朝英，趙樂軍，陳懷琛.紅外TV法眼動測量系統(tǒng)及實現(xiàn)[J].西安電子科技大學學報，1999，26(3):709-712.WANG Chaoying，ZHAO Lejun，CHEN Huaichen.A system and the implementation of eye movement measurement based on the infrared TV[J].Journal of Xidian University，1999，26(3):709-712.

[6]閆國利，田宏杰.眼動記錄技術與方法綜述[J].應用心理學，2004，10(2):55-58.YAN Guoli，TIAN Hongjie.A review of eye movement recording methods and techniques[J].Chinese Journal of Applied Psychology，2004，10(2):55-58.

[7]趙新燦，左洪福，任勇軍.眼動儀與視線跟蹤技術綜述[J].計算機工程與應用，2006，42(12):118-120.ZHAO Xincan，ZUO Hongfu，REN Yongjun.A review of eye tracker and eye tracking techniques[J].Computer Engineering and Application，2006，42(12):118-120.

[8]王軍寧，劉濤，何迪，等.紅外頭盔式眼動議的瞳孔中心定位算法[J].西安電子科技大學學報，2011，38(3):8-12.WANG Junning，LIU Tao，HE Di，et al.Pupil center localization algorithm used for the IR head-mounted eye tracker[J].Journal of Xidian University，2011，38(3):8-12.

[9]周嘉賓.眼動跟蹤系統(tǒng)算法研究與實現(xiàn)[D].西安:西安電子科技大學，2010:241-298.ZHOU Jiabin.Study and implementation on the algorithm of eye tracking system[D].Xi’an:Xidian University，2010:241-298.

[10]閻蓓，王斌，李媛.基于最小二乘法的橢圓擬合改進算法[J].北京航空航天大學學報，2008，34(3):295-298.YAN Bei，WANG Bin，LI Yuan.Optmial ellipse fitting method based on least square principle[J].Journal of Beijing University of Aeronautics and Astronautics，2008，34(3):295-298.

[11]羅軍輝，馮平.MATLAB 7.0在圖像處理中的應用[M].北京:機械工業(yè)出版社，2005:252-257.