基于图像处理技术的瞳孔和角膜反射中心提取算法

王锦榕，袁学海，刘增良

(1.大连理工大学控制科学与工程学院，辽宁大连116024;2.中国人民解放军国防大学信息作战与指挥训练教研部，北京100091)

人总是通过五官从周围世界获取信息，并根据这些信息指导自己的行为，而其中80%～90%的信息是由视觉来获得的，因此，对如何利用眼睛获取视觉信息进行研究是至关重要的.眼动测量一直是国内外有关学者关注与研究的热点课题.近年来，随着传感技术和微电子技术的飞速发展，眼动测量技术得到了很大的发展，测量眼球运动的方法已日趋成熟，在工业控制、机器人学和临床医学上有着广泛的应用前景.目前主要的测量方法有眼电图法[1]、角膜反射法[2]、电磁线圈法[3]、红外光电反射法[4]、红外 TV 法[5]、角膜反射法[6]、瞳孔－角膜跟踪法[7]、Purkinje影像追踪法[7]等.在以上各种方法中，瞳孔－角膜跟踪法测量准确、误差小，且对人无干扰，是目前较为理想的一种眼动测量方法.

文献[8]中采用星射线方法获取瞳孔边界点，但在眼睑遮挡部分瞳孔的情况下，易受到眼睑和睫毛的干扰，很难获取准确的瞳孔边界点.基于最小二乘法和随机原理的椭圆拟合方法[9-10]，通过对48个瞳孔边界特征点进行椭圆拟合，考虑到算法的准确性和实时性，执行的次数一般取100～150之间，该算法的执行时间约为8 ms.为了进一步减小瞳孔椭圆拟合时间，以及在眼睑遮挡部分瞳孔的情况下实现瞳孔中心的准确定位，本文提出了一种基于图像处理技术的瞳孔和角膜反射中心提取算法.

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

在红外光源下，用红外滤光片滤去可见光，通过摄像机实时采集人眼图像视频中的人眼图像，见图1(a).图像中有一个很亮的反射光斑，即角膜反射光斑.瞳孔的形状近似为椭圆，颜色较暗，其外围虹膜颜色较亮，瞳孔与虹膜的边界清晰可见.角膜反射中心和瞳孔中心的提取过程可分为6个步骤:图像预处理、获取图像处理区域、提取角膜反射中心、求取自适应最佳阈值确定瞳孔位置和大小、用梯度法获取瞳孔轮廓特征点、使用椭圆拟合瞳孔轮廓确定瞳孔中心.如图1所示.

图1 图像处理结果Fig.1 The results of image processing

1.1 图像预处理

图像在拍摄和传送过程中会不可避免地受到噪声的干扰，因此先要对图像进行平滑处理以抑制噪声，改善图像质量.常用的平滑处理方法有邻域平均法、中值滤波、高斯滤波等.由于高斯滤波能有效地去除服从正态分布的噪声，且不会被无用的高频信号污染，同时保留了大部分所需要的信号，故采用高斯滤波对人眼图像进行平滑处理.

1.2 获取图像处理区域

在图像处理过程中，主要任务是提取角膜反射中心和瞳孔中心，因此只需要对图像中眼睛所在的区域进行处理即可，从而提高处理速度.由于人眼图像中，瞳孔颜色较暗，即灰度值较低，因此可以采用具有自适应于每帧图像变化的阈值分割法，来获取瞳孔所在的区域.

一般把包含某个目标物体的矩形外框称为该物体所在的矩形区域，它由矩形的左上角坐标和右下角坐标表示.获取图像处理区域的算法如下.

1)取一个较低的阈值t，求出灰度值低于t的面积s和所在的矩形区域R.

2)当s小于某个设定门限阈值时，t加1，再求出 s和 R.

3)重复执行步骤2)，当t达到某个设定值，且s小于设定门限值时，则认为眼睛闭着，退出循环;反之，则继续循环，直到s达到设定门限时求出R，并退出.

将R各边扩大80个像数，使瞳孔和角膜反射光斑在R内，令R为图像处理区域，处理结果见图1(b).

1.3 提取角膜反射中心

由图1(a)可知，角膜反射光斑最大、中心最亮，从中心由里往外亮度逐渐变暗，其中还有一些较小的杂质光斑.当取较大的图像二值化阈值时，角膜反射光斑的面积较小，几乎没有杂质点，但提取的角膜反射中心误差较大.随着二值化阈值的减小，角膜反射光斑的面积逐渐变大，同时引入的杂质点也逐渐增加.因此，要使提取的角膜反射中心误差减小，必须使角膜反射光斑面积足够大，且杂质点很少.

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

1)首先取一个较大的图像二值化阈值，求取粗略的角膜反射中心.

2)以这个角膜反射中心为圆心，r为半径的圆形区域作为图像处理区域，其中r的取值要使得角膜反射光斑在圆形区域内.

3)在图像处理区域内，取一个较小的二值化阈值，求取角膜反射中心坐标(xc，yc).

此方法增大了角膜反射光斑面积的同时，消除了杂质光斑的干扰，既提高了角膜反射中心的提取精度，同时又减少了处理时间，处理结果见图1(c).

1.4 求取自适应阈值确定瞳孔位置和大小

由于瞳孔较暗，灰度值比较小，因此可以采用具有自适应于每帧图像变化的阈值来提取瞳孔.当二值化阈值很小时，只有灰度值小于阈值的点可以通过，则瞳孔面积和瞳孔所在的矩形区域很小，是整个瞳孔的一小部分.随着阈值的增加，瞳孔面积和瞳孔所在的矩形区域面积逐渐加大.由于瞳孔分布均匀，灰度值大小差不多，故瞳孔面积增大的速度比矩形区域快，即瞳孔面积与矩形区域面积的比值逐渐增大.当阈值达到一定值时，即扩大到瞳孔的边缘时，瞳孔面积增加速度逐渐变慢，而矩形区域面积则会继续增加，若继续增加阈值，则瞳孔面积与矩形区域面积的比值逐渐减小.该比值可以表示为

式中:S0表示瞳孔面积，S表示矩形区域面积.随着二值化阈值的逐渐增大，p是一个先增后减的比率函数，此函数极点所对应的阈值就是提取瞳孔的最佳二值化阈值.通过这个最佳阈值就可以求出瞳孔的位置和大小.

由于人的眼睑张合度的不同，会出现眼睫毛遮住瞳孔的现象，造成瞳孔中心与矩形区域中心位置相差很大.在处理过程中，取一个大于瞳孔实际大小的设定值，当瞳孔所在的矩形区域的长度或宽度大于该设定值时，则认为眼睫毛遮住了瞳孔.由于瞳孔中心与矩形区域中心位置相差很大，因此可以通过瞳孔中心到矩形区域各边的距离来重新设置图像处理区域，从而消除眼睫毛的影响，见图2.

在图像处理区域中，求取自适应最佳阈值的具体步骤如下.

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

2)当t＜55时，t加1，求出新的比率 p1和矩形区域R1，R1由矩形的左上角坐标(x11，y11)和右下角坐标(x22，y22)表示，否则退出.

3)当 p－p1＞0.05且 p＞0.6，则此时的 t为最佳阈值，退出;反之，则继续执行.

4)求R1的长度a和宽度b，若满足a＞110或b＞100，则继续执行;反之，令 p=p1，R=R1返回到步骤2)执行.

5)令 R0=R1，R0由矩形的左上角坐标(x1，y1)和右下角坐标(x2，y2)表示.求出瞳孔中心坐标(x，y)，再求瞳孔中心到矩形区域各边的距离，即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为图像处理区域，t减2，求出p和 R，返回到步骤2)执行.

通过求出的自适应最佳阈值确定瞳孔位置和大小，该方法具有很强的鲁棒性，其处理结果见图1(d).

图2 设置处理区域Fig.2 Set the processing area

1.5 用梯度法提取瞳孔轮廓特征点

在人眼图像中，瞳孔和虹膜具有明显的颜色差异，瞳孔的灰度值较低，虹膜的灰度值较高.因此可以通过灰度的变化率和变化方向来衡量瞳孔与虹膜边界的灰度变化，这里分别用梯度的幅值和方向表示.对于某一方向上一定范围内的连续像素点集合，可以用特定的算子求取任意像素点在该方向上的梯度值，梯度值最大的像素点为边界点[11].

通过自适应最佳阈值法确定的瞳孔中心坐标为(xp，yp)，瞳孔在x方向和y方向上的半径分别为rx和ry，这些值具有很高的精度和稳定性.令

因此，瞳孔轮廓特征点一般在以(xp，yp)为圆心，半径为r的圆附近.用梯度法提取瞳孔轮廓特征点的具体步骤如下.

1)以(xp，yp)为圆心，半径为r0的圆环区域作为图像处理区域，见图3，其中r－7≤r0≤r+7.

图3 圆环区域Fig.3 The ring area

2)如图3，在处理区域内，每一组角度θ和半径ri都可以在圆上确定一个相应的像素点pi.其中ri=r－7+i(i=0，1，…，14)，θ为眼图坐标系中，以(xp，yp)为出发点的射线与横轴之间的夹角，0≤θ＜2π.这样可以求出像素点pi在图像中的坐标为(xi，yi)，其中，

3)将像数点 pi(i=1，2，…，14)与圆心相连，在连线上取pi的相邻像数点pi－1，求取pi在θ方向上的梯度值Ti.设像数点pi的灰度值为Li，则

令梯度值最大的像数点为θ方向上瞳孔与虹膜的边界点.

4)取一组不同的θ，执行步骤2)和3)就可以求出一组瞳孔轮廓的特征点.

通过此方法求出的特征点是一组准确的边界点，其处理结果见图1(e).

1.6 椭圆拟合瞳孔轮廓

提取的瞳孔轮廓特征点是一组离散数据，可以通过椭圆拟合来确定瞳孔中心.基于最小二乘法和随机原理拟合椭圆的方法具有很高的精度，但需要一定的执行时间且鲁棒性较差，并且在实时系统中，容易出现干扰数据.由于用梯度法求取的特征点是一组精确的边界点，故可以采用基于几何参数拟合椭圆的方法，椭圆的标准方程为

式中:椭圆中心坐标为(x0，y0)，长半轴为a，短半轴为b，长轴的转角为θ.在人眼图像中，令θ=0，具体的椭圆拟合步骤如下.

1)求取瞳孔轮廓特征点的几何中心(xp，yp)，以(xp，yp)为椭圆中心，a=rx，b=ry拟合椭圆.

2)统计所有特征点与椭圆的距离小于设定值ε的特征点数目n，当n大于设定值时，则拟合成功，执行步骤5);否则继续执行步骤3).

3)以(xp，yp)为椭圆中心，nmax=0，令 a 的初始值为r－7，终止值为 r+7，步长为1，b的初始值为r－7，终止值为r+7，步长为1，进行双重循环.求取每一组a和b对应的n值，当nmax＜n时，令nmax=n，并保存这一组的a和b.

4)当nmax大于设定值时，则拟合成功，保存的a和b就是最佳的长半轴和短半轴，执行步骤5);反之，拟合失败，并退出.

5)除去与椭圆的距离大于设定值ε的特征点，重新求取特征点的几何中心为瞳孔中心.

对150张图像进行瞳孔椭圆拟合，基于最小二乘法和随机原理拟合椭圆方法和基于几何参数拟合椭圆方法的比较结果如表1.从表1中可以看出，该方法具有很高的拟合精度，且执行时间短，其处理结果见图1(f).

表1 2种拟合算法的结果比较Table 1 Comparison of the results of two algorithms

2 实验结果与分析

瞳孔和角膜反射中心的提取算法在以VC++6.0为平台的P4 2.4GHz主频、512MB内存的PC机上实现.用该算法对人眼图像进行处理，得到准确的瞳孔中心和角膜反射中心以及最优的瞳孔椭圆拟合，见图1(f).在眼睑遮挡部分瞳孔的情况下，该算法也能拟合出正确的瞳孔轮廓，见图1(g)和(h).此算法的平均运行时间为15 ms，完全满足实时系统要求.因此可以通过此算法对从摄像机实时采集的人眼图像进行处理，获取瞳孔中心坐标(xp，yp)和角膜反射中心坐标(xc，yc).瞳孔中心相对于角膜反射中心的偏移量称为角膜－瞳孔坐标(xcp，ycp)，即

眼睛注视运动和平滑跟踪运动[7]所对应的角膜－瞳孔位置坐标和实时曲线见图4.图4(a)和(b)为眼睛注视前方某一目标物体所得到的角膜－瞳孔位置坐标和实时曲线;图4(c)和(d)为眼睛围绕着液晶显示屏外框做平滑跟踪运动所得到的角膜－瞳孔位置坐标和实时曲线.实验表明，在眼睑遮住部分瞳孔的情况下新算法实现了瞳孔中心的准确提取，且满足实时性要求.由角膜－瞳孔坐标的实时曲线可知，提取的瞳孔中心和角膜反射中心在实时系统中具有较强的抗干扰能力.

图4 注视运动和平滑跟踪运动Fig.4 The fixation and smooth pursuit

3 结束语

瞳孔中心是确定视线方向最主要的因素，能否提取稳定、准确的瞳孔中心直接影响视线方向的准确性.本算法采用自适应阈值法确定瞳孔的位置，大大减少了处理的数据量，然后用梯度法提取瞳孔轮廓特征点，避免了因眼睑、眼睫毛干扰产生的虚假特征点，最后基于几何参数法的椭圆拟合提取了准确的瞳孔中心，并有效地减少了椭圆拟合时间.后续将通过角膜－瞳孔坐标确定视线方向，进一步研究眼控鼠标.

[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，FUCHS 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法眼动测量系统及实现[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]周嘉宾.眼动跟踪系统算法研究与实现[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].北京:机械工业出版社，2005:252-257.