CN103996020A - 一种头戴式眼动仪检测方法 - Google Patents

Abstract

本发明公开了一种头戴式眼动仪检测方法，包括角膜发射光斑中心提取；提取瞳孔轮廓的特征点；瞳孔轮廓拟合，本发明根据红外光源的抗干扰性计算出眼睛注视点的位置，再经过跟踪算法拟合出闭眼、眨眼时眼睛所在位置，可精确测量眼动注视方位，识别算法精度高，可满足高精度要求下的眼动注视方位的测量。

Description

一种头戴式眼动仪检测方法

技术领域

本发明涉及一种眼睛检测和跟踪的方法，尤其涉及的是一种头戴式眼动仪检测方法。 

背景技术

传统的眼睛检测和跟踪方法主要利用眼睛在外观和形状上与人脸其它部分的差别.眼睛的特点，如黑色的瞳孔、白色的巩膜、圆形的虹膜、眼角、眼睛的形状等，可以用来区分人眼与其他对象，但由于闭眼、眨眼、眼睛大小和位置的可变性、不同的照明条件、人脸的方向等原因，这些眼睛的特点将减小甚至消失。 

发明内容

本发明的目的在于克服现有技术的不足，提供了一种头戴式眼动仪检测方法，利用瞳孔对近红外光的反射现象定位瞳孔。 

本发明是通过以下技术方案实现的，本发明包括以下步骤： 

（1）角膜发射光斑中心提取 

红外光源照射人眼，在角膜表面形成一个亮点，在眼摄像机拍摄到图像中有一个白色区域即角膜反射，在图像平面内，取角膜反射中心作为参考点，根据每帧眼图给出一个自适应变化的阈值，从而对图像进行二值化得到二值图像，所述的自适应变化的阈值根据下面的比率函数得到，该函数的极点所对应的阈值就是自适应后的最佳阈值， 

ratio为比率，Area_最大光斑是最大光斑的面积，Area_总是图像的总面积； 

根据上式得到的阈值对图像进行二值化处理，把大于该阈值的坐标点像素值设为1，小于该阈值的坐标点像素值设为0，Area_最大光斑的面积与只有角膜表面形成的亮点的大小相仿，因此亮点的像素值都为1，图像其它坐标点的像素值都为0，亮点从背景中分离处理， 

对亮点图像区域求质心（i_c,j_c）作为反射中心： 

$i_c=\underset{p_{\mathrm{ij}}\∈{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\Σ}}{p}_{\mathrm{ij}}\·i/\underset{p_{\mathrm{ij}}\∈{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\Σ}}{p}_{\mathrm{ij}}$

$j_c=\underset{p_{\mathrm{ij}}\∈{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\Σ}}{p}_{\mathrm{ij}}\·j/\underset{p_{\mathrm{ij}}\∈{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\Σ}}{p}_{\mathrm{ij}}$

p_ij是（i，j）点坐标的像素灰度值；i为图像的横坐标，j为纵坐标，blob_purk为亮点图像区域； 

（2）提取瞳孔轮廓的特征点 

通过延展星射线法从图像中一个初始点延伸出数条等角度间隔的射线，然后沿着这些射线向外，对射线上的点逐像素求导；当某一点的灰度导数值最大的时候，该点就是是瞳孔轮廓上的一个特征点； 

（3）瞳孔轮廓拟合 

用椭圆方程来表示瞳孔轮廓曲线： 

Ax²+Bxy+Cy²+Dx+Ey+F=0 

其限定条件为： 

B²-4AC<0 

在步骤（2）的特征点中任意提取六个点，将它们代入瞳孔轮廓曲线方程中，利用最小二乘法得出椭圆的系数，拟合出瞳孔轮廓曲线，然后统计与椭圆曲线距离小于δ个像素值的点的数目，并求出这些点占总特征点数目的百分比，当比值大于70%时，此次拟合成功；否则，重新选择六个点，重复上面过程，直到满足判定；当这个过程重复大于50次时，就自动停止计算，返回拟合失败的消息，由前几帧数据来逼近预测当前帧；拟合成功后，利用椭圆曲线方程的系数，就能根据下式求出瞳孔中心坐标： 

$\left\{\begin{array}{c}x=\frac{\mathrm{BE}-2\mathrm{CD}}{4\mathrm{AC}-B^2}\\ y=\frac{\mathrm{BD}-2\mathrm{AE}}{4\mathrm{AC}-B^2}\end{array}\right..$

所述步骤（2）包括以下步骤： 

a、最佳初始点的选择：计算所有特征点与椭圆曲线距离小于单位像素距离δ的特征点坐标在横轴和纵轴的几何平均值，把这个值作为初始点位置坐标； 

b、从初始点等角度辐射出n条射线，并沿着这些射线逐像素求取灰度值的导数， 

O点为所有射线的初始出发点，O、A、B是其中某条射线上相距单位像素距离δ的三点， θ为射线与水平轴的夹角，已知O点坐标和夹角θ以后，推算出A点的坐标值： 

$\left\{\begin{array}{c}{x}_A=x_o+\mathrm{\&delta;}x\cos \mathrm{\&theta;}\\ y_A=y_o+\mathrm{\&delta;}x\cos \mathrm{\&theta;}\end{array}\right.$

逐像素求取O点处灰度导数值： 

$\mathrm{dev}=\frac{I(x_A,y_A)-I(x_o,y_o)}{\mathrm{\&delta;}}$

同理，计算B点坐标，并用之推导A点处灰度导数值，以此类推，得到从初始点发出的每条射线上各个点的灰度导数值，当灰度导数值为最大值时对应的坐标点就把它对应的点记为特征点； 

c、在以步骤b求得的特征点为顶点，以特征点到瞳孔粗略中心连线为中线，左右各间隔a°的两条边所围成的扇形内，第二次延展星射线来提取特征点；以步骤b所得的特征点为出发点，在包括瞳孔粗略中心的2a°的扇形范围内，等角度的发出5*dev／thresh条射线，其中thresh为特征点判定门限值，并在这些射线上再次逐像素求导以提取特征点； 

d、将两次延展星射线法探测到的所有特征点一一做了记录，求出它们的平均位置作为新的瞳孔中心估计值，当新旧两次瞳孔中心点的距离小于或等于1个像素值时，特征点提取成功；否则，从新的瞳孔中心开始，进行新一轮的特征点提取，直到满足新旧两次瞳孔中心点的距离小于5个像素值这一条件，当重复提取特征点的次数大于10次还没有满足条件时，此时认为人眼处于眨眼状态，人为中止运行，进入下一步。 

如闭眼或眨眼，即拟合失败时，用过去两帧的数据逼近下一帧： 

即用第k-1帧和当前第k帧的位置来预测第k+1帧的眼图数据，即用过去两帧去预测下一帧的位置，设置t1=l、t2=2、t3=3， 

分别对应f(t₁)=f(k-1)、f(t₂)=f(k)、f(t₃)=f(k+1) 

$a_0=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)}{N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)}^2}$

$a_1=\frac{N\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)}{N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)}^2}$

把tl和t2，以及f(t₁)和f(t₂)代入上式的a₀和a₁的两个公式，求得： 

$\left\{\begin{array}{c}{a}_0=2\&CenterDot;f(k-1)-f\left(k\right)\\ a_1=f\left(k\right)-f(k-1)\end{array}\right.$

当t3=3时，计算f(t₃)如下所示： 

f(k+1)=f(t₃)=a₀+a₁·t₃=a₀+3a₁

然后得到预测公式： 

f(k+1)=2f(k)-f(k-1)， 

提取出由上一节输出的前两帧眼图数据，然后代入预测公式进行预测即可。 

本发明相比现有技术具有以下优点：本发明根据红外光源的抗干扰性计算出眼睛注视点的位置，再经过跟踪算法拟合出闭眼、眨眼时眼睛所在位置，可精确测量眼动注视方位，识别算法精度高，可满足高精度要求下的眼动注视方位的测量。 

附图说明

图1是本发明的眼动跟踪流程图； 

图2是眼图视频中的一帧； 

图3是角膜发射光斑中心提取流程图； 

图4是逐像素求导的示意图； 

图5是瞳孔椭圆拟合流程图； 

图6是拟合后记录瞳孔轮廓特征图； 

图7是拟合后提取瞳孔中心的示意图。 

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。 

本实施例的识别流程如图1所示，具体包括以下步骤： 

（1）角膜发射光斑中心提取 

红外点光源照射人眼时，在角膜表面会形成一个亮点，因此在眼摄像机拍摄到图像中有一个白色区域即角膜反射，当红外点光源与用户头部保持相对位置不变时，角膜反射在图像中的位置保持不变，所以在图像平面内，通常取角膜反射中心作为参考点。由图2可知，由于红外点光源在图像中形成的角膜反射区域较小，而且在图像眼眶附近范围内其像素灰度值最高，所以可采用一个固定阈值对图像进行二值化得到二值图像，进一步求解该亮点中心可得角膜反射中心的图像坐标。然而，由于用户存在个体差异，以及环境光的变化，采用一个常量阈值进行二值化是不能满足实用要求的。因此，需要根据每帧眼图的具体情况给出一个 自适应变化的阈值。 

当阈值高时，只有亮度比较高的点可以通过。这时的杂质点比较少，但是角膜反射光斑面积也较小，由此提取出的中心坐标误差较大。随着阂值的降低，可以通过的光斑点慢慢增加，从而使得杂质点也有所增加，但是角膜反射光斑的面积同时也变大了，从而使提取中心坐标的误差减小，精度增加。当阈值降低到一定程度后再继续下降时，杂质点还会继续增多，但是角膜反射光斑的面积没有了明显的增加，提取中心坐标的误差将再度增大。的自适应变化的阈值根据下面的比率函数得到，该函数的极点所对应的阈值就是自适应后的最佳阈值， 

ratio为比率，Area_最大光斑是最大光斑的面积，Area_总是图像的总面积；以阈值为自变量的误差函数先减后增，而上式比率函数是一个先增后减的函数，此函数极点所对应的阈值就是自适应后的最佳值。以最佳阈值为门限二值化眼图后，找到面积最大的反射光斑就是角膜反射光斑，如图3所示。 

角膜反射光斑中心不在二值图像中而在暗瞳图像中求取，因为角膜反射中心是角膜反射区域中最亮的部分，越往边界亮度越低，这样可以减小模糊的边界区域对质心结果的影响.方法如下：角膜发射光斑中心在暗瞳图像中求取，对暗瞳图像区域求质心（i_c,j_c）作为反射中心： 

$i_c=\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}\&CenterDot;i/\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}$

$j_c=\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}\&CenterDot;j/\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}$

p_ij是（i，j）点坐标的像素灰度值；i为图像的横坐标，j为纵坐标，blob_purk为亮点图像区域。 

（2）提取瞳孔轮廓的特征点 

红外光照射下的眼图瞳孔轮廓较明显，即该处的灰度值变化较大，延展星射线法正是利用这一点来实现瞳孔轮廓特征点提取的。通过延展星射线法从图像中一个初始点延伸出数条等角度间隔的射线，然后沿着这些射线向外，对射线上的点逐像素求导；当某一点的灰度导数值最大的时候，则该点就是是瞳孔轮廓上的一个特征点； 

a、最佳初始点的选择，初始点的位置越靠近真正的瞳孔中心，迭代的次数就越少，算法就收敛的速度就越快。因为瞳孔是眼图中灰度值最低的部位，根据经验设置一个较低的门限值，本实施例选用80，将灰度值低于这个门限的点统统认为是瞳孔中的点，求这些点的位置坐标的几何平均值，把这个值作为初始点位置坐标，初始点在算法迭代的过程中迅速向真正的瞳孔中心靠拢，当迭代前后两次的射线出发点距离小于某个阈值时，就认为算法收敛，迭代成功； 

b、从初始点等角度辐射出n条射线，并沿着这些射线逐像素求取灰度值的导数， 

如图4所示，O点为所有射线的初始出发点，O、A、B是其中某条射线上相距单位像素距离δ的三点，θ为射线与水平轴的夹角，已知O点坐标和夹角θ以后，推算出A点的坐标值： 

$\left\{\begin{array}{c}{x}_A=x_o+\mathrm{\&delta;}x\cos \mathrm{\&theta;}\\ y_A=y_o+\mathrm{\&delta;}x\cos \mathrm{\&theta;}\end{array}\right.$

逐像素求取O点处灰度导数值： 

$\mathrm{dev}=\frac{I(x_A,y_A)-I(x_o,y_o)}{\mathrm{\&delta;}}$

同理，可以计算B点坐标，并用之推导A点处灰度导数值，以此类推，得到从初始点发出的每条射线上各个点的灰度导数值，当灰度导数值最大的时候，就把它对应的点记为特征点； 

c、在以步骤b求得的特征点为顶点，以特征点到瞳孔粗略中心连线为中线，左右各间隔a°的两条边所围成的扇形内，第二次延展星射线来提取特征点；以步骤b所得的特征点为出发点，在包括瞳孔粗略中心的2a°的扇形范围内，等角度的发出5*dev／thresh条射线，其中thresh为特征点判定门限值，并在这些射线上再次逐像素求导以提取特征点； 

之所以进行第二次星射线延展，是为了能获得更多的特征点，为下一步的瞳孔轮廓拟合提供丰富全面的数据。在上述的扇形区域内发出射线有更大可能穿过瞳孔边缘，灰度值发生突变，增加有效特征点的数日；但是也会有部分光线没有经过瞳孔区域，这样的光线很可能把与眼睑的交点误判为特征点，形成噪声点。由于扇形区域内发出星射线所产生的有效特征点的数目要远远多于新生噪声点的数目，使得这一步显得格外地重要。保障有效特征点数增加的另一个手段是控制发出星射线数目。从初始点发出的射线不再像第一次一样是固定数目了，而是生成(5*dev／thresh)条射线。dev越大，第一次探测到的特征点是有效点的可能性就越高，由它开始第二次探测得到的特征点的有效性也随之提升。 

d、将两次延展星射线法探测到的所有特征点一一做了记录，求出它们的平均位置作为新 的瞳孔中心估计值，当新旧两次瞳孔中心点的距离小于等于1个像素值时，特征点提取成功；否则，从新的瞳孔中心开始，进行新一轮的特征点提取，直到满足新旧两次瞳孔中心点的距离小于5个像素值这一条件，当重复提取特征点的次数大于10次算法还没有收敛时，就有可能是眨眼遮蔽了瞳孔，造成数据中断；于是人为中止算法运行，输入下一帧进行处理。 

（3）瞳孔轮廓拟合 

在上一步特征点提取的过程中，通过迭代会产生一个粗略的瞳孔中心值。这个值是由所有特征点包括有效特征点和噪声点共同生成的，所以存在误差。而我们这一部分要做的就是通过精确的计算排除掉噪声点的影响。瞳孔的轮廓实际上是一个椭圆，并不是正圆形。用椭圆方程来表示瞳孔轮廓曲线： 

Ax²+Bxy+Cy²+Dx+Ey+F=0 

其限定条件为： 

B²-4AC<0 

如图5所示，在步骤（2）的特征点中任意提取六个点，将它们代入瞳孔轮廓曲线方程中，利用最小二乘法得出椭圆的系数，拟合出瞳孔轮廓曲线，然后统计与椭圆曲线距离小于δ个像素值的点的数目，并求出这些点占总特征点数目的百分比，当比值大于70%时，此次拟合成功；否则，重新选择六个点，重复上面过程，直到满足判定；在比较糟糕的眼图中，如眼睑严重地遮挡瞳孔的情形，这个流程可能无法拟合出正确的瞳孔轮廓，陷入无限的循环中；当这个过程重复大于50次时，就自动停止计算，返回拟合失败的消息，由前几帧数据来逼近预测当前帧；拟合成功后，利用椭圆曲线方程的系数，就能根据下式求出瞳孔中心坐标： 

$\left\{\begin{array}{c}x=\frac{\mathrm{BE}-2\mathrm{CD}}{4\mathrm{AC}-B^2}\\ y=\frac{\mathrm{BD}-2\mathrm{AE}}{4\mathrm{AC}-B^2}\end{array}\right..$

图6和图7为拟合结果示意图，图6为记录瞳孔轮廓特征的示意图，图7为提取瞳孔中心的示意图。 

3、眨眼检测 

眼球运动时常出现眨眼情况，将无法进行瞳孔中心定位，因此要对眨眼进行检测并处理，眨眼检测步骤： 

（1）当找不到亮点时，视为眨眼； 

（2）当找不到粗定位瞳孔时，视为眨眼； 

在眨眼时，采用最佳线性逼近预测算法，这样可以充分利用前面图像处理的结果，减少处理计算量，缩短处理时间。 

线性逼近预测算法是根据前几帧中获取的瞳孔中心位置来预测当前帧的瞳孔中心。 

假设原函数f(t)在N个顺序时刻的测量值为f(t)(i=1,2,3...N)，且f(t)可以用如下公式进行最佳线性逼近： 

Y=a₀+a₁t 

t_i时刻测量值与逼近值之间的误差为：Δε_i=f(t_i)-a₀-a₁t_i，对N个点估计的均方误差为： $E\left({\mathrm{\&Delta;\&epsiv;}}_i\right)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{[f\left(t_i\right)-a_0-a_1{t}_i]}^2$

最佳逼近就是使上式取得最小值。 

采用最小二乘法莱计算逼近函数的系数过程如下：E(Δε_i)是关于自变量a₀和a₁的二元函数，欲使E(Δε_i)最小满足： 

$\left\{\begin{array}{c}\frac{E\left({\mathrm{\&Delta;\&epsiv;}}_i\right)}{{\&PartialD;a}_0}=2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}[f\left(t_i\right)-a_0-a_1{t}_i]=0\\ \frac{E\left({\mathrm{\&Delta;\&epsiv;}}_i\right)}{{\&PartialD;a}_1}=-2\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}[f\left(t_i\right)-a_0-a_1{t}_i]\&CenterDot;t_i=0\end{array}\right.$

整理上式得： 

$\left\{\begin{array}{c}{\mathrm{Na}}_0+\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)\&CenterDot;a_0=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)\\ \left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)\&CenterDot;a_0+\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i\right)\&CenterDot;a_0=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)\end{array}\right.$

即： 

$\left[\begin{array}{cc}N& \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\\ \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i& \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]=\left[\begin{array}{c}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)\\ \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)\end{array}\right]$

利用消元法解出方程： 

$a_0=\frac{\left[\begin{array}{cc}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)& \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\\ \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i-f\left(t_i\right)& \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i\end{array}\right]}{\left[\begin{array}{cc}N& \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\\ \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i& \underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i\end{array}\right]}=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)}{N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)}^2}$

同理： 

$a_1=\frac{N\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)}{N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)}^2}.$

所述步骤（3）中，拟合失败时由前几帧数据来逼近预测当前帧的方法为： 

用第k-1帧和当前第k帧的位置来预测第k+1帧的眼图数据，即用过去两帧去预测下一帧的位置，设置t1=l、t2=2、t3=3， 

分别对应f(t₁)=f(k-1)、f(t₂)=f(k)、f(t₃)=f(k+1) 

把tl和t2，以及f(t₁)和f(t₂)代入上述求得的a₀和a₁公式中，得到： 

$\left\{\begin{array}{c}{a}_0=2\&CenterDot;f(k-1)-f\left(k\right)\\ a_1=f\left(k\right)-f(k-1)\end{array}\right.$

当t3=3时，计算f(t₃)如下所示： 

f(k+1)=f(t₃)=a₀+a₁·t₃=a₀+3a₁

然后得到预测公式： 

f(k+1)=2f(k)-f(k-1)。 

提取出由上一节输出的前两帧眼图数据，然后代入预测公式进行预测即可。 

Claims (3)

1.一种头戴式眼动仪检测方法，其特征在于，包括以下步骤：

（1）角膜发射光斑中心提取

红外光源照射人眼，在角膜表面形成一个亮点，在眼摄像机拍摄到图像中有一个白色区域即角膜反射，在图像平面内，取角膜反射中心作为参考点，根据每帧眼图给出一个自适应变化的阈值，从而对图像进行二值化得到二值图像，所述的自适应变化的阈值根据下面的比率函数得到，该函数的极点所对应的阈值就是自适应后的最佳阈值，

ratio为比率，Area_最大光斑是最大光斑的面积，Area_总是图像的总面积；

根据上式得到的阈值对图像进行二值化处理，把大于该阈值的坐标点像素值设为1，小于该阈值的坐标点像素值设为0，Area_最大光斑的面积与只有角膜表面形成的亮点的大小相仿，因此亮点的像素值都为1，图像其它坐标点的像素值都为0，亮点从背景中分离处理，

对亮点图像区域求质心（i_c,j_c）作为反射中心：

$i_c=\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}\&CenterDot;i/\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}$

$j_c=\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}\&CenterDot;j/\underset{p_{\mathrm{ij}}\&Element;{\mathrm{blob}}_{\mathrm{purk}}}{\mathrm{\&Sigma;}}{p}_{\mathrm{ij}}$

p_ij是（i，j）点坐标的像素灰度值；i为图像的横坐标，j为纵坐标，blob_purk为亮点图像区域；

（2）提取瞳孔轮廓的特征点

通过延展星射线法从图像中一个初始点延伸出数条等角度间隔的射线，然后沿着这些射线向外，对射线上的点逐像素求导；当某一点的灰度导数值最大的时候，该点就是是瞳孔轮廓上的一个特征点；

（3）瞳孔轮廓拟合

用椭圆方程来表示瞳孔轮廓曲线：

Ax²+Bxy+Cy²+Dx+Ey+F=0

其限定条件为：

B²-4AC<0

在步骤（2）的特征点中任意提取六个点，将它们代入瞳孔轮廓曲线方程中，利用最小二乘法得出椭圆的系数，拟合出瞳孔轮廓曲线，然后统计与椭圆曲线距离小于δ个像素值的点的数目，并求出这些点占总特征点数目的百分比，当比值大于70%时，此次拟合成功；否则，重新选择六个点，重复上面过程，直到满足判定；当这个过程重复大于50次时，就自动停止计算，返回拟合失败的消息，由前几帧数据来逼近预测当前帧；拟合成功后，利用椭圆曲线方程的系数，就能根据下式求出瞳孔中心坐标：

$\left\{\begin{array}{c}x=\frac{\mathrm{BE}-2\mathrm{CD}}{4\mathrm{AC}-B^2}\\ y=\frac{\mathrm{BD}-2\mathrm{AE}}{4\mathrm{AC}-B^2}\end{array}\right..$

2.根据权利要求1所述的一种头戴式眼动仪检测方法，其特征在于，所述步骤（2）包括以下步骤：

a、最佳初始点的选择：计算所有特征点与椭圆曲线距离小于单位像素距离δ的特征点坐标在横轴和纵轴的几何平均值，把这个值作为初始点位置坐标；

b、从初始点等角度辐射出n条射线，并沿着这些射线逐像素求取灰度值的导数，

O点为所有射线的初始出发点，O、A、B是其中某条射线上相距单位像素距离δ的三点，θ为射线与水平轴的夹角，已知O点坐标和夹角θ以后，推算出A点的坐标值：

$\left\{\begin{array}{c}{x}_A=x_o+\mathrm{\&delta;}x\cos \mathrm{\&theta;}\\ y_A=y_o+\mathrm{\&delta;}x\cos \mathrm{\&theta;}\end{array}\right.$

逐像素求取O点处灰度导数值：

$\mathrm{dev}=\frac{I(x_A,y_A)-I(x_o,y_o)}{\mathrm{\&delta;}}$

同理，计算B点坐标，并用之推导A点处灰度导数值，以此类推，得到从初始点发出的每条射线上各个点的灰度导数值，当灰度导数值为最大值时对应的坐标点就把它对应的点记为特征点；

c、在以步骤b求得的特征点为顶点，以特征点到瞳孔粗略中心连线为中线，左右各间隔a°的两条边所围成的扇形内，第二次延展星射线来提取特征点；以步骤b所得的特征点为出发点，在包括瞳孔粗略中心的2a°的扇形范围内，等角度的发出5*dev／thresh条射线，其中thresh为特征点判定门限值，并在这些射线上再次逐像素求导以提取特征点；

d、将两次延展星射线法探测到的所有特征点一一做了记录，求出它们的平均位置作为新的瞳孔中心估计值，当新旧两次瞳孔中心点的距离小于或等于1个像素值时，特征点提取成功；否则，从新的瞳孔中心开始，进行新一轮的特征点提取，直到满足新旧两次瞳孔中心点的距离小于5个像素值这一条件，当重复提取特征点的次数大于10次还没有满足条件时，此时认为人眼处于眨眼状态，人为中止运行，进入下一步。

3.根据权利要求2所述的一种头戴式眼动仪检测方法，其特征在于，拟合失败时，用过去两帧的数据逼近下一帧：

即用第k-1帧和当前第k帧的位置来预测第k+1帧的眼图数据，即用过去两帧去预测下一帧的位置，设置t1=l、t2=2、t3=3，

分别对应f(t₁)=f(k-1)、f(t₂)=f(k)、f(t₃)=f(k+1)

$a_0=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)}{N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)}^2}$

$a_1=\frac{N\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;f\left(t_i\right)-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\&CenterDot;\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}f\left(t_i\right)}{N\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t^2}_i-{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{t}_i\right)}^2}$

把tl和t2，以及f(t₁)和f(t₂)代入上式的a₀和a₁的两个公式，求得：

$\left\{\begin{array}{c}{a}_0=2\&CenterDot;f(k-1)-f\left(k\right)\\ a_1=f\left(k\right)-f(k-1)\end{array}\right.$

当t3=3时，计算f(t₃)如下所示：

f(k+1)=f(t₃)=a₀+a₁·t₃=a₀+3a₁

然后得到预测公式：

f(k+1)=2f(k)-f(k-1)，

提取出由上一节输出的前两帧眼图数据，然后代入预测公式进行预测即可。