基于眼震影像的多方向上眼震信号提取方法

摘要

本发明涉及一种基于眼震影像的水平性眼震、垂直性眼震和旋转性眼震信号的提取和分析方法，包含的步骤是：解码和归一化；检测和补偿反光点；虹膜区域粗定位；计算亮度梯度分布矩阵；构建梭形三叉树；计算最优路径；检测瞳孔边界和虹膜边界；提取水平及垂直眼震信号；提取虹膜特征并编码及提取旋转眼震信号。本发明以眼震影像为处理对象，通过视频分析和数字图像处理技术，提取其眼震信号，并以曲线图的方式直观地展现在计算机屏幕上，提供医学上有益的参考信息。

说明书

所属技术领域

本发明涉及视频影像分析及数字图像处理领域，特别涉及一种以眼震影像为对象，对其进行视频解码、图像分割、特征提取、特征编码和模板匹配，获取眼震视频中的水平性眼震、垂直性眼震和旋转性眼震的提取和分析方法。

背景技术

眼震，又称眼球震颤(Nystagmus)，是一种不受意念控制的眼球节律性运动，前庭性眼震由交替出现的慢相和快相运动组成。眼震症状与中枢神经系统疾病有着密切的关联，例如老年痴呆、癫痫、脑肿瘤、前庭周围病变、中枢性病变(包括脑干、小脑病变)、眩晕症和某些眼病。因此，通过观察眼震可以了解以上疾病的病变情况，甚至可以作为诊断眩晕症和前庭病变的主要依据。眼震按方向分为水平性眼震、垂直性眼震和旋转性眼震。

目前对于眼震的检查，除了裸眼观测外，最主要技术为眼震电图检测法，通过在眼球两侧放置电极，采集眼球运动产生的生物电信号，通过信号放大器的作用放大该信号，再以图形的方式表示出来，这就是眼震电图(ElectroNystagmoGraphy,ENG)。眼震电图检测法可定量分析眼睁、闭眼时的眼震频率和眼震振幅，在临床检查中起到了重要的作用。但是在试验中，不同人的角膜视网膜电位常常不同，甚至同一个人在不同时间或不同情绪的情况下也表现不同，所以每次检测时都需要定标且容易受病人情绪的干扰。另外由于表示眼震的生物电信号极其微弱，放大器大约需要放大2万倍才能检测出明显的信号，其他生物信号的干扰或者环境中的电磁信号污染均可对检测结果造成影响。另一个局限性是目前该技术只能检测水平性眼震和垂直性眼震，无法检测旋转性眼震。而旋转性眼震是眩晕症和前庭病变的主要临床病征。此外，眼震电图检测对于设备和环境要求很高，且设备价格高昂、不易普及，目前国内仅在大型医院中使用。

眼震视图检测法(VideoNystagmoGraphy,VNG)由于检测方便、成本低廉且不受到病人主观情绪因素和环境电磁信号的干扰，逐渐受到国内外医学机构的重视。VNG系统采用红外线照明和高速摄像头进行眼震影像的采集，通过观察和分析眼震影像来获得水平性眼震、垂直性眼震和旋转性眼震。但是对于VNG的相应研究还非常少，临床上还是基本上依靠肉眼观察眼震影像来进行诊断，未能有可靠的自动化眼震分析和提取技术得以普及，主要由以下因素造成：(1)影像复杂。采集过程中的外界因素如光照强度、摄像头位置和个体本身的行为特征如眼球转动、眼睑覆盖、眨眼、闭眼等均反应在采集到眼震影像中。而眼震电图中没有存在如此复杂的信号；(2)成像畸变。摄像头通过光学系统将物体投影到平面感光元件上从而获得影像。而人体的眼球是由角膜、巩膜、虹膜、晶状体、玻璃体等组织构成的复杂类球形的器官。采集过程中，三维的眼球信息投影到二维的感光元件上造成信息的截取和丢失。常见的问题有：A.眼球的旋转引起的眼球光学主轴与摄像头光学主轴不共线形成投影畸变；B.摄像头光学主轴不过眼球的几何中心形成的离轴畸变。处理影像过程中必须采用适应性较好的模型才能避免以上畸变对检测的影响。(3)研究很少。目前国内外对于眼震影像的分析还处于起步阶段，研究成果基本上是应用现有虹膜识别系统中的虹膜分割和提取的技术，而由于在采集过程中虹膜识别系统对人眼的行为具有高度约束，且采集的过程持续时间很短。故能避免影像出现(1，2)中所提出的干扰因素，但是眼震视图检测中则无法避免。现有的虹膜分割技术采用的道格曼(Daugman)微积分算子、霍夫变换(Hough)、主动轮廓模型(ActiveContour)等方法均在复杂的眼球影像上表现较差，造成系统鲁棒性较差，特别是有畸变发生时，故目前未获大量应用，主要还是依赖医师主观判断。

发明内容

本发明的目的是针对现有眼震视图检测法中应用的眼震分析和提取方法的不足，结合机器学习、图像分割、特征提取和模式识别技术，提出了一种基于眼震影像的多方向上眼震信号提取方法。该方法可以从复杂眼震影像中分割瞳孔和虹膜，并进行特征提取、编码和匹配，提取三个方向上的眼震信息，并以眼震曲线图的方式展现在计算机屏幕上，提供医学上有益的参考信息。

本发明涉及包含以下步骤：

a)解码和归一化：在红外LED光源照射下，将已经采集好的眼震影像文件或者实时采集的眼震影像视频流按照其编码方式解码为图像帧，取待处理的帧，将其转化为灰度模式，并保持其长宽比s，放大或压缩原图像至大小为(W×W·s)的归一化图像，W限定归一化后的图像分辨率，根据计算机性能设定。

b)检测和补偿反光点：遍历上述归一化后的眼球图像的每个像素，若像素的亮度值I大于反光点候选亮度阈值T_I,则认为该像素具有反光点的亮度特征，然后对通过亮度阈值判定的像素检测其圆形邻域的亮度梯度特征判断其是否为反光点像素。检测方法为统计检测半径r从r_min递增至r_max的过程中，计算以像素为中心、半径为r的圆形邻域内所有像素的亮度的累加值，寻找累加值对半径r的偏导数的最小值γ，若γ小于反光点候选梯度阈值T_G则判定该像素在反光点区域，否则说明该像素为亮度较高的人体部位，例如巩膜。其中[r_min,r_min]为反光点可能的半径范围,检测过程的计算公式如下：

其中G_σ为高斯滤波器，x₀,y₀为待测像素坐标

γ(x₀,y₀)实际反应了以该像素为中心的邻域扩大过程中，邻域区域的亮度最大递减梯度，即表征邻域向外扩张时亮度的衰减剧烈程度，而反光点边界正是存在剧烈的亮度衰减。通过反光点候选梯度阈值T_G筛选出拥有该特征的像素，则该像素既满足亮度特征又满足梯度特征，那么认为该像素为反光点像素。然后按照坐标对所有检测出的反光点像素进行聚集和分类，将邻近的像素合并为区域，从而获得若干反光点区域进行填补。

由于红外LED光源是一个点光源，其镜面反射的成像是一个圆形点，故反光点区域近似一个高亮的小圆。填补过程首先计算出反光点区域的几何中心(x_r,y_r)和以该中心为圆心可以扫过区域内所有反光点像素的最小半径r，构成圆然后在圆的外围做环圆C为能够包裹该区域的最小圆，环R为圆C的外接环作为信息域为填补圆C内反光点像素提供亮度信息，环R内的像素称为信息点，见图2。根据邻近相似理论，与待填补点距离越近的信息点对其产生的影响越大，反之越小。故以待填补点与圆心连线在环R内所经过的像素的亮度值为变量，像素与圆心的距离所占以上所有所经过像素与圆心距离之和的比例为概率，所求得的期望值为该待填补点的亮度值。简化的计算方式如下：

$I(x,y)=\frac{D_1(x-x_r)}{{2\mathrm{rr}}^{\prime }}\underset{i=1}{\overset{r}{\Sigma }}\left(s(i,\mathrm{ki})I(x_r+i,F(y_r+\mathrm{ki}))\right)+\frac{D_2(x-x_r)}{{2\mathrm{rr}}^{\prime }}\underset{i=-1}{\overset{-r}{\Sigma }}\left(S(i,\mathrm{ki})I(x_r+i,F(y_r+\mathrm{ki}))\right),\left|k\right|\le 1$

$I(x,y)=\frac{D_1(x-x_r)}{{2\mathrm{rr}}^{\prime }}\underset{i=1}{\overset{r}{\Sigma }}\left(s(\mathrm{ki},i)I(F(x_r+\mathrm{ki}),y_r+i)\right)+\frac{D_2(x-x_r)}{{2\mathrm{rr}}^{\prime }}\underset{i=-1}{\overset{-r}{\Sigma }}\left(S(\mathrm{ki},i)I(F(x_0+\mathrm{ki}),y_0+i)\right),\left|k\right|>1$

$S(x,y)=\left(\begin{array}{c}1,\mathrm{ifr}<\sqrt{x^2+y^2}\le r+r^{\prime }\\ 0,\mathrm{else}\end{array}\right)$

$D_1\left(x\right)=\left(\begin{array}{c}{d}_{\mathrm{near}},x\ge 0\\ d_{\mathrm{far}},x<0\end{array}\right),D_2\left(x\right)=\left(\begin{array}{c}{d}_{\mathrm{far}},x\ge 0\\ d_{\mathrm{near}},x<0\end{array}\right)$

其中k为待填补点与圆心连线的斜率，d_near为连线上待填补点到环的较近距离，d_far为较远距离，F为向下取整函数，函数S确保只有在信息域内的像素才能提供亮度信息。显然信息域内的所有像素都将被用到补偿反光点过程中。

c)虹膜区域粗定位：对经过步骤b)消除反光点后的图像输入到预先训练好的AdaBoost分类器中，定位虹膜区域。如果有区域通过分类器判定，则从原始图像裁剪出该区域，命名为关键区域，并计算关键区域的中心坐标(x_k,y_k)。反之，若无任何区域通过判定，则认为该图像内不存在明显的瞳孔和虹膜，是由于采集过程中闭眼或眨眼造成，则输出闭眼或眨眼标记后结束本帧处理，跳转到步骤a)处理下一帧。

d)计算亮度梯度分布矩阵：以坐标(x_k,y_k)为中心将关键区域图像由直角坐标系转化到极坐标系(ρ,θ)下，其中ρ为极径，θ为极角，获得极坐标图像。构造梯度检测算子P用于分析极坐标图像的亮度梯度分布，梯度检测算子是一个长度为n像素的线性单调递增信号，斜率通过分析大量的虹膜图像获得，取值为统计虹膜图像的瞳孔边界和虹膜边界向外扩张时亮度递增梯度的平均值。

亮度梯度的检测方法为，遍历极坐标图像的每一个像素，以其坐标为中心，在极径ρ方向上计算信号P与其邻域的亮度乘积的累加值作为该像素的检测结果，从而获得梯度分布矩阵。由于在瞳孔边界和虹膜边界处均存在类似的亮度递增梯度，该梯度检测算子可以放大该边界处的亮度递增信号，削弱其他区域的亮度不变信号和递减信号，故在梯度分布矩阵中瞳孔边界和虹膜边界所对应位置上将出现局部峰值，而此局部峰值将成为检测边界的主要依据。梯度分布矩阵中的每个元素对应极坐标图像上的每个像素，其元素值则表示对应像素极径方向上的亮度梯度。

e)构建梭形三叉树：由眼球图像的特性分析可以得出，瞳孔边界和虹膜边界沿极坐标图像极角θ方向生长，可以抽象为f(θ,ρ)和h(θ,ρ)的连续光滑曲线。由于在关键区域限定的区域内，瞳孔边界区域和虹膜边界区域含有最强的亮度递增梯度，故在反应梯度变化的梯度分布矩阵中，边界曲线为矩阵中所经过元素的累加值最大的连续光滑路径。检测该路径的方法为，以极角θ递增方向为矩阵的行方向，极径ρ递增方向为列方向，以矩阵的第一列的元素为根节点，沿行方向，构建梭形三叉树，从而可以构建出ρ个树，见图4。每个三叉树包含从其根节点出发满足连续光滑性质的所有路径。生成的三叉树集合为{Tree_i|0≤i≤ρ-1}。

f)计算最优路径：探测三叉树从根节点到叶子节点所经过的节点值之和最大的路径，认为该路径为该树的最优路径，记录为R_i＝{(x₀,y₀),(x₁,y₁)…(x_θ-1,y_θ-1)}，其节点值之和为V_i。故通过计算每一个三叉树的最优路径，可以获最优路径和节点累加值的序列集合{(R_i,V_i)|0≤i≤ρ-1}。每个序列代表以第i行第一列的元素为起始点，所经过的元素累加值最大的连续光滑路径和及其累加值。

g)检测瞳孔边界和虹膜边界：遍历序列集合{(R_i,V_i),0≤i≤ρ-1}，寻找累加值V最大的序列(R_pupil,V_pupil),由e)中的结论，路径R_pupil为极坐标下的瞳孔边界曲线f(θ,ρ)离散化后的像素集合。从序列集合中移除起始点位于R_pupil以下的所有序列，再次寻得累加值最大的序列(R_iris,V_iris)，则路径_Riris为极坐标下的虹膜边界曲线h(θ,ρ)离散化后的像素集合。由于梭形三叉树模型其紧凑的节点间关系，在该模型下能保证边界的光滑且检测过程可塑性强，更能够适应复杂畸变的情况,见图4。

h)提取水平及垂直眼震信号：将上一步中获得的路径R_pupil所有像素映射到直角坐标系下，获得关键区域内虹膜的轮廓边界点集合R′_pupil并计算该轮廓的几何中心(x_pupil,y_pupil)。假设当前处理的帧为眼震影像的第t帧，则有目前眼球的水平位置信息和垂直位置信息则水平位移矢量为${\mathrm{\Delta x}}_{\mathrm{pupil}}^t=x_{\mathrm{pupil}}^t-x_{\mathrm{pupil}}^{t-1}$及垂直位移矢量为${\mathrm{\Delta y}}_{\mathrm{pupil}}^t=y_{\mathrm{pupil}}^t-y_{\mathrm{pupil}}^{t-1}.$在眼震曲线图中，水平眼震信号和垂直眼震信号分别对该组位移矢量，追加该时刻的位移矢量到图中并绘制曲线，得到从0到t帧的水平眼震信号曲线l_t和垂直眼震信号曲线v_t。

i)提取虹膜特征并编码：在极坐标图中虹膜边界R_pupil和瞳孔边界R_iris之间的图像信息为虹膜带。截取虹膜带，并归一化为360×100大小，归一化过程使确保虹膜信息可以均匀排列并使极角θ的精确度为1°，见图5。然后通过Gabor滤波器对其进行特征提取，并对其进行编码，从而获得该帧的虹膜码矩阵M_t。

j)提取旋转眼震信号：以上一帧的虹膜码矩阵M_t-1为窗口，通过在水平滑动窗口，检测与当前虹膜码矩阵M_t重叠区域相似度最大时的水平滑动位移Δr。相似度的衡量依据为重叠区域内元素的海明距离的平均值，该值越小表示相似度越高。由于通过归一化后，极角θ的精确度为1°，眼球不可能在短短的一帧时间内旋转超过20°，故约束窗口滑动的最大位移为20。在眼震曲线图中，旋转眼震信号对应该旋转位移矢量Δr，该值为正代表逆时针旋转，负为顺时针，图中追加该时刻的旋转位移矢量并绘制曲线，得到0到t帧的旋转眼震信号曲线r_t，见图6。跳转到步骤a)处理后续图像帧。

本发明专利的有益效果是：(1)对眼震影像的处理及分析以获得水平眼震信号、垂直眼震信号和旋转眼震信号，比常规方法更完备；(2)能处理眨眼、闭眼、反光、畸变等复杂干扰因素，稳定性和适应性更强；(3)采用独创的梭形三叉树模型，比传统模型检测更为精确、计算量更小。

附图说明

图1是本发明专利的眼震影像提取眼震信号的处理流程图

图2反光点填补环形模型结构

图3梭形三叉树模型结构

图4瞳孔边界和虹膜边界检测结果

图5归一化后的虹膜带

图6眼震信号曲线图

具体实施方式

下面结合一个非限定的实施例对本发明作进一步说明

参见处理流程图

本发明的算法采用C语言结合开源的数字图像处理开发库OpenCV在Windows平台上实现，图形界面和曲线图绘制由C#语言在.NET4.0框架下实现。界面功能包括眼震影像的展示和回放、同步的眼震信号曲线图、影像路径选择、全局参数设定等。

主要模块的内容介绍如下：

(1)对眼震影像进行解码和归一化：由于眼震视图检测在国内未广泛普及，实验采用的眼震影像由实验设备采集。采集设备介绍如下：由于在红外照射下眼球会出现亮瞳效应，本发明实验时采用8个波长为850nm的红外LED组成环形光源实现照明。同时采用SONY-HD960H高清CCD板机和微距镜头进行眼震影像采集。照明模块和采集模块均放置于完全不透光的眼罩中，通过佩戴就可实现舒适的检测。摄像头采集的速度为25帧/秒，采用MPEG-2编码，帧分辨率为720×576。取当前需要处理的帧，通过MPEG-2解码器进行解码后，转化为灰度图并压缩为分辨率640×512的归一化图像。

(2)检测和补偿角膜反光点：为了保持对不同光源亮度的适应和算法处理速度，试验中设定反光点候选亮度阈值T_I为220，寻找图像中亮度大于T_I的像素进一步分析。取检测半径范围为[r_min＝2,r_max＝5]依据分析公式对其进行圆形邻域的亮度梯度特征分析，获得通过判定的反光点像素。我们采用k-means算法对获得的反光像素进行聚类，获得若干反光点区域。然后对每个反光点区域进行填补，设定环的宽度为5构建环形信息域，按照b)中的填补公式完成反光点补偿。

(3)虹膜区域粗定位：本发明中采用1190张手动分割的虹膜区域图片为正样本，5000张不包含眼球的图片为负样本进行AdaBoost分类器训练。每一张正样本图片都归一化为64×64大小，并提取局部二值模式(LocalBinaryPatterns)特征，计算公式如下：

$\mathrm{LBP}(x_c,y_c)=\underset{n=0}{\overset{7}{\Sigma }}s(i_n-i_c)2^n$

其中(x_c,y_c)为待提取特征的像素又称为中心点，遍历该中心点的8个相邻像素，i_c和i_n分别代表中心点的亮度和相邻像素的亮度，s则是一个符号函数。试验中训练输出的分类器包括19层和50个特征。应用该预先训练的分类器在反光点填补后的图像中查找是否能通过该19层判定的窗口，如果存在则截取该窗口图像命名为关键区域，否则说明不包含明显的瞳孔区域，跳转至(1)处理下一帧。

(4)计算亮度梯度分布矩阵：以上一步获得的关键区域几何中心为中心，转化为极坐标图像。取长度为12斜率为0.5构造梯度检测算子P，试验中取值为[-3,-3,-2,-2,-1,0,0,1,2,2,3,3]，按照d)中计算规则应用于极坐标图像中获得亮度梯度分布矩阵。

(5)检测瞳孔边界和虹膜边界：假设亮度梯度矩阵行数为Num_rows则依照梭形三叉树构建规则将获得三叉树集合为{Tree_i,0≤i≤Num_rows-1}。计算每颗树的最优路径和最大元素累加值，获得序列集{(R_i,V_i),0≤i≤Num_rows-1}。在序列集中寻找累加值V最大的序列(R_pupil,V_pupil)，则路径R_pupil为极坐标下瞳孔边界曲线f(θ,ρ)离散化后的像素集合。在移除路径起始点位于R_pupil下的所有序列，获得截取后的序列集{(R_i,V_i),0≤i≤pupil-1}，寻找V值最大序列(R_iris,V_iris)，则路径R_iris为极坐标下的虹膜边界曲线h(θ,ρ)离散化后的像素集合。

(6)绘制水平及垂直眼震信号：将上一步中获得的路径R_pupil所有像素映射到直角坐标系下，则直角坐标系下瞳孔边界像素点的集合为{(x_i,y_i)|0≤i<n},并计算该边界的几何中心：

$x_{\mathrm{pupil}}=\frac{{\Sigma }_{i=0}^n{x}_i}{n},y_{\mathrm{pupil}}=\frac{{\Sigma }_{i=0}^n{y}_i}{n}$

取上一帧进行比较，计算水平位移矢量及垂直位移矢量并绘制在眼震曲线图中，见图6。

(7)绘制旋转眼震信号：极坐标图中虹膜边界曲线h(θ,ρ)和瞳孔边界曲线f(θ,ρ)之间的图像信息为虹膜带。截取虹膜带，压缩并归一化为360×100大小。采用Gabor滤波器对归一化的虹膜带进行特征提取，Gabor滤波器的构造方法和试验中的参数取值如下：

实数部分：$g(x,y;\lambda ,\theta ,\psi ,\sigma ,\gamma )=e^{(-\frac{x^{\prime 2}+{\gamma }^2{y}^{\prime 2}}{{2\sigma }^2})}\cos (2\pi \frac{x^{\prime }}{\lambda }+\psi )$

虚数部分：$g(x,y;\lambda ,\theta ,\psi ,\sigma ,\gamma )=e^{(-\frac{x^{\prime 2}+{\gamma }^2{y}^{\prime 2}}{{2\sigma }^2})}\sin (2\pi \frac{x^{\prime }}{\lambda }+\psi )$

其中x′＝xcosθ+ysinθ,y′＝-xsinθ+ycosθ，取θ值为0滤波方向为水平滤波，σ取值为0.5，λ为Ψ为0，γ为1。，然后根据滤波后的每个像素的实部值和虚部值的符号对其进行编码，具体编码规则为：

从而获得该帧的虹膜码矩阵M_t。取上一帧的虹膜码M_t-1为窗口，约束窗口滑动的最大位移为20，在位移区间[-20,20]内滑动窗口，计算每次位移窗口M_t-1与M_t重叠区域内元素的海明距离的平均值。取海明距离的平均值最小时的滑动位移Δr，则其为相似度最大时的水平滑动位移，将Δr绘制在眼震曲线图中，见图6。如果所有帧处理完毕，分析过程结束并保存眼震曲线图，否则跳转至(1)步骤处理一下一帧。

以上结合流程图所描述的具体实施例仅是本发明在所采用的实验环境中的优选实施方式，并非对本发明专利保护范围的限定，任何基于本发明精神所做的改进都在本发明保护范围之内。