刘 帅,刘元宁,朱晓冬,冯家凯 ,卢时旭
1.吉林大学 符号计算与知识工程教育部重点实验室,长春 130012
2.吉林大学 软件学院,长春 130012
3.吉林大学 计算机科学与技术学院,长春 130012
虹膜因其具有的唯一性、非侵犯性、稳定性以及天然防伪性[1]等特点,使得准确率较其他生物特征高。因此,虹膜识别技术是当前比较热门的生物特征识别技术之一。
虹膜识别过程分为虹膜图像采集,虹膜质量评价,虹膜定位,虹膜归一增强,虹膜特征表达与匹配[2]。虹膜定位就是定位出虹膜内外边界,将两个边界间的虹膜环形区域分割出来,以方便虹膜特征表达和识别[3]。传统的虹膜定位算法有Daugman提出的微积分圆模板法[4]。基于边缘检测的hough法[5]。冯薪桦博士运用二值图像下边缘粗定位虹膜内边界,并利用内边界圆心左右两侧的一维灰度均值信号粗定位虹膜外边界[6]。Tan提出采用边缘检测和圆拟合的算法实现虹膜定位[7]。Yu采用基于小范围搜索的虹膜定位算法[8]。这些算法虽然能准确定位虹膜的位置,但鲁棒性较差,对图像要求都较高,并且对眼毛眼睑,光照等噪声干扰的抗干扰能力较差,识别效果不稳定。
因此,本文提出基于分块搜索的虹膜定位算法。由内到外、由粗到细对虹膜进行定位。用多种虹膜库在相同识别算法下进行处理,实验证明,该算法可以有效增强虹膜定位的鲁棒性,并且不同采集器在不同条件下所采集的虹膜图像均可以使用该算法进行虹膜定位,识别效果稳定。
2.1.1内圆去光斑处理
虹膜内圆就是眼睛的瞳孔部分,大多数虹膜采集器出于增强光照以及定位需要,会在摄像头旁加装LED灯。因此实际采集的眼睛图像中瞳孔会有一圈光斑,光斑会给虹膜定位造成极大的干扰,因此首先要去除光斑的影响。
首先根据眼睛中瞳孔、眼毛、虹膜、巩膜、皮肤之间灰度值的差异,用阈值处理法[9]将眼睛图像转为二值图像。眼睛图像中,瞳孔和眼毛灰度值最小,虹膜灰度值比它们大,因此,选取灰度图中第一个波峰之后的第一个波谷的值T作为阈值[10]。图1的灰度直方图如图2所示,阈值化公式如式(1)所示:
图1 采集的眼睛图像
图2 眼睛灰度直方图以及阈值选取点
B(x,y)代表阈值化后的灰度值,H代表点(x,y)的实际灰度值。图1的图像阈值化结果如图3所示。
图3 眼睛图像的二值化图像
由图3可知,虽然瞳孔与部分眼毛已经被二值化为白色,但瞳孔内部还是有四个黑色的LED灯光斑。为去掉这些光斑。先用Sobel算子将瞳孔内的光斑边缘检测出来,Sobel算子是图像处理中用于边缘检测的一种常见算子[11]。之后对光斑图像依次进行开运算和闭运算[12]。运算公式如式(2)所示:
圆形结构元素用符号b和c表示。眼睛图像用符号 f表示。得到的结果图像如图4所示。
图4 光斑轮廓图像
将得到光斑轮廓与二值化图像叠加,并将叠加后的图像用种子填充算法[13]填充瞳孔内光斑,进而得到无光斑二值化图像。如图5所示。
图5 无光斑眼睛二值化图像
2.1.2 内圆粗定位
用没有光斑的二值化图像做虹膜内圆定位,首先用基于边缘检测的hough法[14]进行内圆粗定位。hough圆检测就是利用图像的特征检测出目标轮廓,将检测到的连续或不连续的点连接起来,呈现出一种几何形状。利用这一特点,对无光斑二值化图像进行圆检测。hough变换过程图公式(3)所示:
表示圆的参数空间。
当g(xj,yj,xc,yc,r)=0,表示以(xc,yc)为圆心,以r为半径的圆通过了边界点(xj,yj)。当H出现最大值的时候,表示边界点经过该圆的时候,数量最多,认定该圆为边界的轮廓圆。
本文对无光斑的二值化图像进行处理,根据具体虹膜库设定的参数寻找一个合适大小的圆认定为瞳孔所在圆。因为二值图像中还是有眼睑眼毛等干扰,可能造成圆的位置不准,偏离瞳孔。因此需要将所圈的圆的位置的周围点的灰度值与灰度直方图中的瞳孔灰度值比较,进而判断所圈圆的位置是否为瞳孔。图1的边缘检测图如图6所示。内圆粗定位图如图7所示。圈定截图如图8(a)所示。圈定图灰度直方图如图8(b)所示。
图6 无光斑眼睛二值化图像边缘检测图
图7 虹膜粗定位图
图8 圈定区域截图及其灰度直方图
因为hough变换中参与计算的参数较多,而每张虹膜图像的采集环境是不同的,因此同样的参数不可能满足所有虹膜,虽然hough圆检测法可以找到瞳孔大致位置,但仍然无法完全匹配上瞳孔的轮廓。
2.1.3运用分块搜索法的内圆精定位
精定位之前,首先要判断粗定位定位出的瞳孔范围与正式的瞳孔范围相比是偏大还是偏小,还是粗定位的内圆已经达到精确定位的要求。
采用分块搜索法,先根据瞳孔粗定位时得到的圆心和半径,在无光斑眼睛二值化图像四个方向上边界邻域上分别划分正方形边界检测算子(检测算子的大小根据具体虹膜库而定,图1虹膜的检测算子为20×20)。检测图如图9所示。
图9 最初的分块检测算子
设定一个偏小阈值和一个偏大阈值,分别计算四个检测算子范围内的平均灰度值(分别为Nleft,Nright,Nup,Nunder),如果四个值中有的值过小,小于偏小阈值,说明这个方向上的瞳孔圈偏小,需要扩大瞳孔范围。反之,如果四个值中有的值过大,大于偏大阈值,说明圈的范围过大,需要缩小。本例中Nleft,Nup,Nunder的值全部为0,Nright为142.8,说明这上、下、左三处在瞳孔内部,圈的范围偏小,需要扩大圈定瞳孔范围。
之后在四个方向边界邻域上设置移动算子(移动算子的大小根据具体虹膜库而定,图1虹膜的移动算子竖直方向为10×20,水平方向为10×20)。本例最初状态如图10所示。
图10 四边界移动算子最初状态
根据需要平移的方向,不断平移移动算子,每移动n个像素点,就重新划分新的检测算子,检测算子大小不变,并计算检测算子的平均灰度值d,当d的值大于偏小阈值的时候(如果圈的范围过大就看是否小于偏大阈值),记下移动距离m,横纵坐标向该方向平移m的一半,半径缩小(或放大)m的一半。
设粗定位时划定的圆心为(xc,yc),半径为r。先统计需要平移的方向的个数。若只需要向一个方向平移。以向左移动为例,当移动算子在该方向上平移距离为m时,圆心坐标为,半径为;若需要向两个方向平移,以圆心坐标分别两个方向平移,半径根据三角函数得出,以左上平移为例,左平移距离为m,上平移距离为n时,圆心坐标,半径为如果需要向三个方向平移,则先按某一方向调整圆心坐标与半径,调整完毕后,再沿另外两个方向平移,以左上下方向为例,移动算子先向左平移m,得到新圆心坐标为,新半径为,之后再同时调整上下移动算子,上算子移动k,下算子移动t。最终的圆心坐标为最终半径为若四个方向均需要调整,先在左右方向平移调整移动算子做平移m,右平移n,新圆心坐标为新半径为 r±,之后再同时调整上下移动算子,上算子移动k,下算子移动t。最终的圆心坐标为最终半径为
图1图像中左上下三个方向均需要调整,调整示意图如图11所示。
图11 左侧方向调整示意图
图1 的精定位内圆图像如图12所示。
图12 内圆精定位
2.2.1 外圆粗定位
虽然虹膜的内外边界不是同心圆,但具有一定的耦合关系[15]。因此外圆粗定位是结合内圆精定位得到的瞳孔圆心进行外圆粗定位。
首先,对图像使用一个5×3的算子进行卷积操作,对卷积后的结果进行二值化操作。卷积操作图像如图13(a)所示,二值化图像如图13(b)所示。
图13 卷积操作及其二值化
计算规定角度范围弧线上的白点的数量,需要针对不同的半径都要计算,半径的变化范围从(r,2×r)开始逐渐变大。由于虹膜图像上下部分容易受到眼睑眼毛等噪音的干扰,为了减小噪音对虹膜外边界定位的干扰,因此仅计算一定虹膜角度下的白点数量,由图13(b)可以看出,虹膜右下侧白点较多,易读取。为了在受干扰较小的情况下尽可能多地读取白点的数量,选取角度范围为瞳孔圆心水平方向正右侧顺时针0°~60°[16],之后按公式(5)计算:
图14 外圆粗定位
2.2.2运用分块搜索法的外圆精定位
卷积操作的外圆粗定位存在两个问题,一是定位的外圆是以瞳孔为圆心的同心圆。二是因为每张图像拍摄的条件不同,导致无法准确预估虹膜半径与瞳孔半径的比值。只能根据生理学的实践得到的人眼虹膜半径约为瞳孔半径1~2倍这一结论[17],将最大半径设为瞳孔半径的2倍。因此,需要再次使用分块搜索法进行外圆精确定位。
外圆的分块与内圆类似,但亦有不同,需要根据外圆粗定位的圆心和半径在水平方向两边界处划分检测算子(大小根据具体虹膜库改变,本文例子使用10×20),如图15所示。
图15 外圆粗定位检测算子示意图
根据检测算子的灰度直方图,左右移动检测算子,将检测算子灰度直方图中的灰度分布与眼睛图像灰度直方图中的灰度分布进行比较,分别取左右两侧灰度值分布达到巩膜区域内的灰度值(第二波谷左右),检测算子向左移动距离k,向右移动距离为m。粗定位圆心为(xc,yc),半径为R。移动后,圆心为半径为精确定位的外圆如图16所示。
图16 外圆精定位
本文使用了吉林大学自主采集第三代虹膜库,选取了其中的60类别,每类别选取20张图像,共计1 200张图像;中科院CASIA-Iris-interval虹膜库,选取了其中50类别,每类别选取了15张,共计750张图像。图像事先经过质量评价检测[18]合格,但其中包含了多种类型,包括贴近边界,较重眼毛等类型图像,用于检测算法鲁棒性。
实验将本文的虹膜定位算法与Daugman的微积分圆模板法[4]、hough圆定位法[5]、文献[8]的小范围搜索法进行比较。
本文中,四种方法统一使用Daugman教授提出的橡皮圈模型[19]进行虹膜图像归一化,将环形虹展开成一个512×64的矩形大小。并截取其中纹理最强的部分,截成256×32固定矩形区域。之后增强[20]图像纹理。图1的全环矩形如图17(a)所示。截取图像如图17(b)所示。截取增强图像如图17(c)所示。
图17 图像归一,增强示意图
本实验使用的CPU主频为双核2.5 GHz,内存为8 GB,操作系统为Windows XP sp3。在特征提取前,先对所有归一化虹膜图像进行水平移位消除虹膜旋转问题[21]。之后,统一采用Daugman博士提出的二维Gabor滤波进行处理[22],筛选并且提取特征纹理信息,通过比较样本间的欧式距离[23]判定虹膜图像的归属。
ROC曲线[24]是一种表示错误拒绝率(False Reject Rate,FRR)以及错误接收率(False Accept Rate,FAR)关系的曲线。ROC曲线反映了虹膜识别算法的匹配性能。其越接近横纵坐标轴,就代表系统的性能就越好。当FRR与FAR的值相等的时候,就是系统性能最好的时候,这个值称为等错率(Equal Error Rate,EER)。EER越小,性能越好。另外,正确识别率(Correct Recognition Rate,CRR)[25]也是常用的虹膜识别系统的性能评定指标之一。
本文根据CRR、EER、ROC曲线以及算法在内圆平均时间(Pupil Time,PT);外圆定位平均时间(Iris Time,IT)和整体平均定位时间(T),以及定位算法成功定位的虹膜数量进行算法的比较。
吉大虹膜库实验结果如表1所示,ROC曲线如图18所示。CASIA-Iris-interval实验结果如表2所示,ROC曲线如图19所示。
表1 吉林大学虹膜库实验结果
表2 CASIA-Iris-interval实验结果
图18 吉林大学虹膜库实验结果的ROC曲线
图19 CASIA-Iris-interval虹膜库的实验结果的ROC曲线
本文提出了一种基于分块搜索的虹膜定位算法,先用阈值处理法将眼睛图像转化为二值图像,然后利用sobel算子消除图像中定位光斑的影响,用基于边缘检测的hough算法粗定位内圆;用分块搜索法精定位内圆;用卷积操作粗略定位外圆;再利用分块搜索精确定位虹膜外圆。使用两种不同的虹膜库,与另三种定位方法在相同特征表达以及识别算法下进行比较。实验结果表明,相较于传统定位算法,本文方法的ROC曲线更接近于横纵坐标轴,平均耗时有所缩短并且准确率有所提升,抗干扰性有所加强,具有良好鲁棒性。