CN101093538A - 一种基于小波变换过零表示的虹膜识别方法 - Google Patents

Abstract

本发明公开了一种基于小波变换过零表示的虹膜身份识别算法，此算法能够在欠理想的情况下，即所采集的虹膜图像存在一定眼皮遮挡和一定形变的，仍然能够通过小波变换过零表示的方法准确的提取特征点，并进行虹膜编码，从而完成正确的身份识别。同时，本发明对原来图像进行了基于一定阀值的二值化处理，所以能够提高瞳孔粗定位的精度和简化瞳孔精定位，并能提高精度。提出基于虹膜外缘粗定位的双阀值处理，进行中值滤波后，能够有效提取虹膜的大概部分，并进行快速比较准确的粗定位。在之后的Hough变换精定位过程中要使用粗定位的结果，并且比较精确的定位虹膜外缘。使用4阶B－样条小波滤波器进行特征提前并进行虹膜编码，有效的提取信息部分，简化编码，提高匹配的效果。

Description

一种基于小波变换过零表示的虹膜识别方法

技术领域

本发明属于生物识别验证技术领域，特别涉及生物识别验证技术领域中的虹膜识别验证技术领域。

背景技术

在高度信息化的现代社会，人们的身份鉴别已经渗透到日常生活的每一个方面。同时，由于交通、通讯和网络技术的飞速发展，人类的活动范围越来越大，身份鉴别的难度和重要性也越来越突出。对于我们这样一个人口众多的国家，身份鉴别有着特别广泛的应用前景和重要的战略意义，金融、安全、网络、电子商务等无一不需要可靠的身份鉴别。传统的利用密码等身份鉴别的方法具有易遗忘、易假冒等缺点，已不符合现代数字社会的需求。基于生物特征的身份鉴别技术利用人本身所拥有的生物特征来判别人的身份，这些生物特征具有“人各有异、终生不变、随身携带”三个特点，具有稳定、便捷、不易伪造等优点，近年来已成为身份鉴别的热点。常用的生物特征包括指纹、掌纹、虹膜、脸像、声音、签名和笔迹等。

虹膜识别技术是最近发展起来的身份鉴别技术，具有识别准确性高、速度快，具有防伪性和非侵犯性等优点。虹膜识别已有一定规模的商业应用，如机场检票系统、自动柜员机(ATM)等，另外在网上支付、门禁保安、远程登录以及对罪犯或嫌疑人的识别与管理等领域都有很好的应用前景。这种包括虹膜识别在内的基于人体生物特征的身份鉴别技术已经形成了一个专门的学科，称为生物特征识别。

1936年，眼科专家Frank Burch指出虹膜具有独特的信息，可用于身份识别。1987年，眼科专家Aran Safir和Leonard Flom首次提出了利用虹膜图像进行自动虹膜识别的概念，但是他们并没有开发出一个实际的应用系统。1991年，美国洛斯阿拉莫斯国家实验室的Johnson发明了一个虹膜识别系统，这是有文献记载的最早的一个应用系统。1993年，John Daugman发明了一个高性能的虹膜识别原型系统，目前，大部分自动虹膜识别系统使用Daugman的核心识别算法。1996年，Richard Wildes研制成功基于虹膜的身份认证系统。但是这些成功研制的虹膜识别系统都要求有精密的光学前端系统输入装置，以保证质量良好的、没有严重上下眼皮遮挡的、没有形变的虹膜图像。但他们所使用的虹膜库都是非流通的，目前所常见的免费可使用的库只有中科院的CASIA1.0、2.0虹膜库。但这个库中的大部分虹膜图像存在上下眼皮遮挡问题，部分存在弹性形变。

现在通常使用的虹膜特征方法有：

(1)Daugman提出的基于Gabor滤波器的虹膜纹理相位编码算法，该算法利用Gabor滤波器的局部性和方向性对虹膜纹理进行分解编码，要求输入图像质量高。详见文献：J.G.Daugman，How Iris Recogniton Works，IEEE Trans on Circuits & Systems for VideoTechnology，volume 14，issue 1，2004，21-30.

(2)Wildes提出用高斯型滤波器在不同分辨率下分解虹膜图像，并把结果进行存储比较。上述两种算法的性能容易受外界环境影响，眼球表面反光、虹膜的几何变形、睫毛遮挡、噪声等问题也严重影响识别性能，且实际应用中对测试环境和设备要求较高。详见文献：Wildes，R.P，Iris recognition：an emerging biometric technology，Proceedings of the IEEE，Volume85，Issue 9，Sept.1997，Page(s)：1348-1363。

(3)Boles提出的基于小波变换的虹膜识别方法：它有效地克服了图像平移、旋转、图像缩放等前者未能解决的问题，提高了识别性能，但具体的方法描述不清楚。详见文献：W.W.Boles，B.Boashash，A human Identification Technique Using Images of the Iris and WaveletTransform.IEEE Trans，on signal processing，volume 46，issue 4，1998，1185-1188.

以上3种方法的实验都是通过他们自己的虹膜库来验证效果的，在不同的其他实验条件下很难实现，而且识别效果最好的方法(1)Daugman方法是建立在高质量图像(很小的眼皮遮挡和形变、高清晰度)的基础上的，在现在没有标准虹膜库的条件下，很难完全使用这些方法。

发明内容

本发明提供了一种基于小波变换过零表示的虹膜身份识别算法，此算法能够在欠理想的情况下，即所采集的虹膜图像存在一定眼皮遮挡和一定形变的，仍然能够通过小波变换过零表示的方法准确的提取特征点，并进行虹膜编码，从而完成正确的身份识别。流程图如说明书附图1所示。

为了方便描述本发明地内容，首先将相关概念说明如下：

1.灰度图像、灰度值：灰度图像是指图像中只包含了亮度信息而没有任何颜色信息的图像。本专利所使用的灰度图像是8级(2⁸)灰度，即取值是0～255。一般将白色的灰度值定义为255，黑色灰度值定义为0，而由黑到白之间的明暗度均匀地划分为256个等级。

2.瞳孔和虹膜外缘的表示：因为它们都近似于圆，所以可以分别用圆心和半径来表示它们的位置和大小。瞳孔圆心表示为(x₀，y₀)、半径表示为r；虹膜圆心表示为(x_i，y_i)、半径表示为R。

3.二值化处理：设定的阈值T0，对于原来图像I中点(x，y)的灰度值I(x，y)，小于T0的取为0，大于T0的取为1。说明书附图2为二值化阈值选取的示意图。说明书附图3(a)为采集的虹膜图像，说明书附图3(b)为二值化后的虹膜图像效果图。

4.迭代运算，通过r，x₀，y₀在一定范围内的取值，这里是x₀＝x_p+α，α＝-10～+10；y₀＝y_p+β，β＝-10～+10；r＝r_p+k*Δr，k＝-10～+10，其中瞳孔粗定位结果圆心(x_p，y_p)、半径r_p。计算后面的表达式所得到最大(最小)值时的x₀，y₀，r的取值。

5.中值滤波：所谓中值滤波，是指把以某点(x，y)为中心的小窗口内的所有象素的灰度按从大到小的顺序排列，将中间值作为(x，y)处的灰度值(若窗口中有偶数个象素，则取两个中间值的平均)。

6.Canny算子边缘检测：Canny边缘检测法利用高斯函数的一阶微分，它能在噪声抑制和边缘检测之间取得较好的平衡。具体检测步骤见文献(《MATLAB6.X图像处理》，清华大学出版社，247页～250页)：

(1)用高斯滤波器来对图像滤波，可以去除图像中的噪声。

(2)用高斯算子一阶微分对图像进行滤波，得到每个像素梯度的大小|G|和方向θ。

$\left|G\right|={[{\left(\frac{\∂f}{\∂x}\right)}^2+{\left(\frac{\∂f}{\∂y}\right)}^2]}^{\frac{1}{2}}$

$\mathrm{\θ}={\tan }^{-1}[\frac{\∂f}{\∂y}/\frac{\∂f}{\∂x}]$

f为滤波后的图像。

(3)对梯度进行“非极大抑制”。

梯度的方向可以被定义为属于4个区之一，各个区用不同的邻域像素用来进行比较，以决定局部极大值。这4个区及其相应的比较方向如下表3-1所示。

表格3-1四个区域其相应的比较方向

4	3	2
			1	X	1
2	3	4

例如，如果中心像素x的梯度方向属于第4区，则把x的梯度值与它的左上和右下相邻像素的梯度值比较，看x的梯度值是否是局部极大值。如果不是，就把像素x的灰度设为0，这个过程称为“非极大抑制”。

(4)对梯度取两次阈值得到两个阈值T1和T2，T1＝0.4*T2。我们把梯度值小于T1的像素的灰度设为0，得到图像1。然后把梯度值小于T2的像素的灰度设为0，得到图像2。由于图像2的阈值较高，去除了大部分噪声，但同时也损失了有用的边缘信息。而图像1的阈值较低，保留了较多的信息。我们可以以图像2为基础，以图像1为补充来连接图像边缘。

(5)连接边缘的具体步骤如下：

①图像2进行扫描，当遇到一个非零灰度值的像素P时，跟踪以P为开始点的轮廓线，直到该轮廓线的终点Q。

②考察图像1中与图像2中Q点位置对应的点Q’的8-邻域区域。如果Q’点的8-邻域区域中有非零像素R’存在，则将其包括到图像2中，作为点R。从R开始重复①步，直到我们在图像1和图像2中都无法继续为止。

③当完成对包含P的轮廓线的连接之后，将这条轮廓线标记为已访问。回到第①步，寻找下一条轮廓线。重复步骤①、②、③，直到图像2中找不到新轮廓为止。

7.Hough变换：Hough变换是图像处理中边缘连接的一种常用的方法。

虹膜具有显而易见非常良好的环状特性，其边缘为圆形边缘。x，y空间平面中的每个圆都可以用其圆心坐标(x₀，y₀)和半径r三个参数表示为：(x-x₀)²+(y-y₀)²＝r²，

根据Hough变换的基本思路，x，y空间平面中的每个圆对应x₀，y₀，r参数空间上的一个点。对x，y空间平面中的每个边缘点(x_i，y_i)，我们可以给所有与该点对应的x₀，y₀，r空间的直方图方格一个增量。当所有的边缘点实行完这种操作后，包含(x₀，y₀，r)的方格将具有局部最大值。然后对x₀，y₀，r空间的直方图进行局部最大值搜索可以获得边界圆的参数。

8.光照不均的校正：先取16×16大小图像中平均值作为图像背景的照度；然后将粗略估计出的背景照度矩阵扩展成和原始图像大小相同的矩阵，这可以通过双三次插值实现；将估计出的背景照度从原始图像中减去，即可修正照度不均的影响，但这样做的后果是图像变暗；可以通过调整图像的灰度来进行校正，最后得到进行了光照不均的校正后的图像。

9.双三次插值：是一种更加复杂的插值方式，它能创造出比双线性插值更平滑的图像边缘。通过双三次插值创造了一个象素，而这个象素的象素值是由它附近的(4×4)个邻近象素值推算出来的，因此精确度较高

10.直方图均衡化：所谓直方图就是在某一灰度级的象素个数占整幅图像的象素比h＝nj/N，其中nj是灰度级在j的象素数，N是总象素数，扫描整幅图像得出的h的离散序列就是图像的直方图，h求和必然等于1，所以直方图可以看成是象素对于灰度的概率分布函数。

直方图是高低不齐的，因为象素灰度是随机变化的，直方图均衡化就是用一定的算法使直方图大致平和。算法如下：

对于一个直方图，设Pr(r)是原始图像直方图，Ps(s)是均衡化的直方图，由于其是一个概率分布函数，所以有Ps(s)ds＝Pr(r)dr，因为要进行均衡化，令Ps(s)＝1，得ds＝Pr(r)dr/1两边积分得s＝∫Pr(r)dr数字图像是离散的，因此离散化上式得 $\mathrm{sk}=\underset{j=0}{\overset{k}{\mathrm{\Σ}}}\mathrm{nj}/N$ 左式k，j是离散量下标，N是象素总数

由此得出每一象素的sk为均衡化后的正规化灰度(即灰度正规化到[0，1])，统计sk即可得出均衡化后的直方图。

在均衡化过程中可以对每一象素映射到新的实际灰度值sk*255，就实现了图像的变换。

11.虹膜图像归一化：虹膜归一化的目的是将每幅原始图像调整到相同尺寸和对应位置，从而消除平移、缩放和旋转对于虹膜识别的影响。假设已经得到虹膜的内、外边界圆周的参数分别为：(x₀，y₀，r₀)、(x_i，y_i，R)，我们采用极坐标变换方法进行归一化。

由于虹膜内外圆边界不是同心的，所以这种极坐标变换也不是同心的。

设内圆圆心为I(I_x，I_y)，半径为R_i；外圆圆心为O(O_x，O_y)，半径为R₀。根据I(I_x，I_y)，O(O_x，O_y)的相对位置，可以分为几种情况：1)I_x＝O_x；2)I_x＞O_x，I_y≥O_y；3)I_x＞O_x，I_y≤O_y；4)I_x＜O_x，I_y≥O_y；5)I_x＜O_x，I_y≤O_y。不失一般性，以2)为例子，如说明书附图5所示，以内圆圆心作为极坐标系统的中心，做与水平线成H角的射线，它与内、外边界各有一个交点，分别记做B(x_i(θ)，y_i(θ))，A(x₀(θ)，y₀(θ))。则有： $\mathrm{\α}=\mathrm{arctg}\left|\frac{I_y-O_y}{I_x-O_x}\right|,$ ∠OIA＝π-θ+α， $\∠\mathrm{OAI}=\arcsin \frac{R_i\·\sin \∠\mathrm{OIA}}{R_0},$ ∠IOA＝π-∠OIA-∠OAI， $\mathrm{IA}\left(\mathrm{\θ}\right)=\sqrt{{R_i}^2+{R_0}^2-2\·R_i\·R_0\·\cos \∠\mathrm{IOA}}$ 其它四种情况的结果可以类似推导出来。

于是，射线上两个交点之间的任何一点都可以用A(x₀(θ)，y₀(θ))，B(x_i(θ)，y_i(θ))的线性组合表示：

x(r，θ)＝(1-r)*xi(θ)+r*x0(θ)，

y(r，θ)＝(1-r)*yi(θ)+r*y0(θ)

其中r∈[0，1]，θ∈[0，2π]。

该变换将虹膜图像中的每个点一一映射到极坐标对(r，θ)中去.这种由直角坐标下的虹膜图像到极坐标下的映射可以表示为：

I(x(r，θ)，y(r，θ))→I(r，θ)

这种映射对于平移和内外圆环的大小变换具有不变性。而对于旋转变化，将在后面的算法中提取旋转不变的特征。综上所述，整个虹膜识别的过程对于虹膜图像的旋转、平移和尺度变化以及瞳孔的收缩都具有不变性。归一化之后的图像在极坐标系(r，θ)下，展开成一个大小为64×512矩形，如说明书附图4(f)。

一种基于小波变换过零表示的虹膜身份识别算法，其特征在于包括如下步骤，如图1所示：

步骤1，将虹膜摄像采集装置靠近人眼，在人眼进入摄像范围内后进行图像采集，通过数字图像采集CCD的光电转换得到含有虹膜图像的原始灰度图像I(x，y)，图像左上角为I(x，y)的零点I(0，0)，图像向右令横坐标x递增，图像向下令纵坐标y递增，得到如说明书附图3(a)所示的人眼图像，虹膜图像的原始灰度图像I(x，y)的取值范围为0～255，x的范围为1～M，y的范围为1～N，原图像大小为M×N；

步骤2，显示原始灰度图像I(x，y)的直方图，如说明书附图2所示，标识1左面的第一个峰值为瞳孔区域的灰度值，取阈值T0为其向右的第一个谷值，如图3中1所示。

步骤3，对步骤1中得到的原始灰度图像I(x，y)，用阈值T0进行二值化处理。

原始灰度图像I(x，y)小于T0的点I(x，y)的值取为0，I(x，y)大于T0的点I(x，y)的值取为1，组成二值矩阵I’(x，y)，得到的二值化图像I’(x，y)大小等于原图像I(x，y)，如说明书附图3(b)所示。

步骤4，瞳孔粗定位。

通过步骤3得到原始灰度图像I(x，y)的二值化后的图像I’(x，y)，用式子 $P\left(x\right)=\underset{y}{\mathrm{\Σ}}{I}^{\′}(x,y),$ $P\left(y\right)=\underset{x}{\mathrm{\Σ}}{I}^{\′}(x,y)$ 将二值化后的图像进行水平方向上和垂直方向上的投影，得到二值化后的图像水平方向上和垂直方向上的投影值P(x)，P(y)；并且通过计算 $X_p=\arg \underset{x}{\min }P\left(x\right),Y_p=\arg \underset{y}{\min }P\left(y\right),$ 得到瞳孔粗定位圆的圆心在直角坐标上的结果(x_p，y_p)，如说明书附图3(c)所示。从水平的投影值P(x)谷值点X_p，向左搜索，得到点P1，满足在P1以后的点的投影值变化比较平缓(在5以内)，同理，从谷值点X_p向右搜索得到P2，满足在P2以后的点的投影值变化比较平缓(在5以内)。取P1、P2之间部分的范围为瞳孔的粗定位直径，并且可以得到瞳孔粗定位圆的半径r_p＝P2-P1。

步骤5，瞳孔精定位方法：根据步骤4的瞳孔粗定位结果，圆心(x_p，y_p)、半径r_p，通过迭代运算求得瞳孔精确定位结果，圆心(x₀，y₀)、半径r，方法为：

$\max (r,x_0,y_0)|\underset{m}{\mathrm{\Σ}}{I}^{\′}(x_{m,k},y_{m,k})-\underset{m}{\mathrm{\Σ}}{I}^{\′}(x_{m,k-1},y_{m,k-1})|,---\left(1\right)$

其中X_m，k＝r×cos(m×Δθ)+x₀，x₀＝x_p+α，α＝-10～+10，

Y_m，k＝r×sin(m×Δθ)+y₀，y₀＝y_p+β，β＝-10～+10，

r＝r_p+k×Δr，k＝-10～+10，α、β、k为变量

Δθ是圆周上的一个角度采样间隔，取Δθ＝1°，所以m＝0～360，Δr为在半径方向上的间隔，取Δr＝1。

通过α、β、k的代入运算，求得当式子(1)得到最大值时的x₀、y₀、r的值，则瞳孔的精确定位结果：圆心为(x₀，y₀)、半径为r₀。

步骤6，虹膜外缘粗定位方法：从步骤5得到的瞳孔精定位结果，由虹膜外缘半径范围的先验知识，取一个与瞳孔同心、半径值为比瞳孔半径值大20个像素的同心圆，求出这个圆和瞳孔圆之间的虹膜部分的灰度值平均值Avg：

Avg＝位于这环形区域点的灰度值的和÷位于这环形区域之间的总点数；

计算过程是：

环形区域中的点的集合s：{(x，y)|瞳孔圆心²≤x²+y²≤(瞳孔圆心+20)²}原图像上的灰度值为I(x，y)，位于这环形区域点的灰度值的和为 位于这环形区域之间的总点数为集合s的点的总数。

令阈值T1＝Avg-15，阈值T2＝Avg+15，作为双阈值对原灰度图像I(x，y)进行如下处理：

这样得到图像矩阵I”(x，y)是基本上只包含虹膜部分的图像，再对其进行中值滤波，得到图像矩阵I₂(x，y)。

求I₂(x，y)在水平投影值 $P^{\′}\left(x\right)=\underset{y}{\mathrm{\Σ}}{I}_2(x,y),$ 然后求 ${X_p}^{\′}=\arg \underset{x}{\min }{P}^{\′}\left(x\right),$ 从谷值点X_p’向左搜索，得到点P1’，满足在P1’以后的点的点间灰度值变化比较平缓(在5以内)；同理，从谷值点X_p’向右搜索得到P2’，满足在P2’以后的点的点间灰度值变化比较平缓(在5以内)。取P1’、P2’之间部分的范围为虹膜外缘的粗定位直径，并且可以得到虹膜外缘粗定位圆的半径R’＝P2’-P1’，虹膜外缘粗定位圆心设为与瞳孔精定位圆心一样，为(x₀，y₀)。

步骤7，虹膜外缘精定位方法：对步骤3二值化后的图像，进行基于Canny算子的边缘提取，得到关于边缘的二值图像，1为边缘，0为其他点。根据步骤6的虹膜外缘粗定位半径R’和圆心(x₀，y₀)，在R’±δ，x₀±α，y₀±β范围内和在图像的-45°～45°and 135°～225°范围内进行Hough变换，其中δ、α、β为待定常数，其取值与精度和运算量有关，它们的值越大精度越高、运算量越大，反之亦然。这样就可以得到虹膜外缘的精定位结果：半径为R，圆心为(x_i，y_i)。

步骤8，对原图像进行归一化处理：由步骤5的瞳孔精定位结果(x₀，y₀，r₀)和步骤7的虹膜外缘精定位结果(x_i，y_i，R)，进行虹膜图像归一化处理。把原灰度图像的虹膜部分归一化为64×512的矩形图像I_nor(i，j)，1≤i≤64、1≤j≤512，如说明书附图3(f)。

(详见J.G.Daugman，“High confidence visual recognition of persons by a text ofstatistical independence”.IEEE Trans.on Pattern Analysis and Machine Intelligence，Volume15，issue 11，1993，1148-1161.和王蕴红、朱勇、谭铁牛，《基于虹膜识别的身份鉴别》，自动化学报，Vol.28，No.1，Jan，2002)

步骤9，对步骤8得到的图像进行光照不均的校正，然后进行直方图均衡化，得到总体亮度均匀、纹理清晰的图像I_adj(i，j)，1≤i≤64、1≤j≤512，如说明书附图3(h)。

步骤10，虹膜编码：对步骤9所得到的图像进行特征点提取。先用4阶B-样条小波滤波器对步骤9所得到的图像I_adj(i，j)分别在水平方向和垂直方向进行卷积滤波。4阶B-样条小波滤波器：是一个小波滤波器，其参数如下：

表1 4阶B-样条小波滤波器

K是变量，ψ(K)是相应的小波函数值。卷积方法如下：

temp_h(i，j)＝ψ(-i)*I_adj(i，j)……行方向上卷积，j固定

temp_v(i，j)＝Ψ(-j)*I_adj(i，j)……列方向上卷积，i固定

分别得到行和列方向上的卷积滤波后的临时存储矩阵temp_h(i，j)，temp_v(i，j)，并且进行如下判决得到行和列方向上编码，行和列方向上的编码进行与操作得到虹膜编码：

3.令code(i，j)＝code_h(i，j)And code_v(i，j)，这里的And为求两个数的与运算，以上i＝1～64，j＝1～512。

这样二值化矩阵code(i，j)就是由虹膜特征点所组成的虹膜编码，大小为64×512。编码情况的例子如说明书附图5所示。

步骤11，编码匹配：由步骤9得到的矩阵code(i，j)可以代表原始图像的虹膜的编码，如果要把两个虹膜编码code1(i，j)和code2(i，j)进行匹配，可以使用下面的Hamming距离的式子，其表征的意义为两个虹膜编码之间的相似程度：

$\mathrm{Ham}\min g\_\mathrm{Dis}\tan \mathrm{ce}=\underset{\mathrm{\φ}}{\mathrm{Min}}\left\{\frac{1}{M\×N}\underset{i=0}{\overset{M}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}\left(\mathrm{code}1(i+\mathrm{\φ},j)\right)\mathrm{Xor}\left(\mathrm{code}2(i,j)\right)\right\}$

for-10≤φ≤10

如果HammingDistance小于某阈值T，则两个虹膜是相同虹膜；如果HammingDistance大于T，则两个虹膜是不相同虹膜。其中T的取值和样本的性质有关，可以在数据统计的时候决定最优的阈值T，如说明书附图6所示。

通过上面步骤，我们可以得到每一个采集到图像的虹膜编码，并且可以通过匹配来分析是同一虹膜还是不同虹膜。在实际应用中，只要在数据库中记录虹膜编码就行了，在身份识别的时候只要将采集到的虹膜进行出来并编码后，与在库中的虹膜编码进行匹配就行了。

需要说明的是：

步骤1所采集的原始图像一般为灰度图像，如果是彩色图像则需要通过灰度转换把彩色图像转换为灰度图像，具体方法是：对于彩色RGB模型，灰度值＝0.3×红色+0.59×绿色+0.11×蓝色。

通过以上步骤，本发明得的有益效果是：

1.对原来图像进行了基于一定阈值的二值化处理，所以能够提高瞳孔粗定位的精度和简化瞳孔精定位，并能提高精度。

2.提出基于虹膜外缘粗定位的双阈值处理，进行中值滤波后，能够有效提取虹膜的大概部分，并进行快速比较准确的粗定位。在之后的Hough变换精定位过程中要使用粗定位的结果，并且比较精确的定位虹膜外缘。

3.使用4阶B-样条小波滤波器进行特征提前并进行虹膜编码，有效的提取信息部分，简化编码，提高匹配的效果。

附图说明

图1是总的流程图；

图2为二值化阈值选取的示意图，1表示选取的阈值T0；

图3：(a)是CASIA 1.0的一幅原始图像；(b)分离瞳孔的二值化图像；(c)瞳孔二值化图像水平和垂直灰度值投影；(d)在阈值T1和T2之间的虹膜分离图；(e)精定位的效果图；(f)归一化后的虹膜图像；(g)归一化虹膜图像进行光照估计补偿后的图像；

图4为归一化的示意图，2表示虹膜外边界的圆心为O(O_x，O_y)，3表示虹膜内边界的圆心为I(I_x，I_y)，4表示外边界上的点A(x_o，y_o)，5表示内边界上的点B(x_i，y_i)，6表示外边界半径为R_o，7表示内边界半径为R_i，；

图5：(a)虹膜特征空间的分布图(虹膜1)(b)虹膜特征空间的分布图(虹膜2)(白色点为虹膜特征点)；

图6分别为相同类之间和不同类之间的虹膜的Hamming距离分布图，8表示相同类之间Hamming距离分布，9表示不同类之间Hamming距离分布；

图7采集的各种虹膜图像；

图8为虹膜瞳孔精定位效果图。

具体实施方式

采用本发明的方法，在机场和门禁系统入口处前，和其他任何合适位置采用虹膜摄像装置采集认证人的虹膜图像；其次把采集到的虹膜原始图像做为原数据输入到用Matlab语言编写的软件中进行处理；最终得到一个认证人的虹膜识别编码。

实验共采用108*7＝756张CASIC 1.0虹膜库图像做为原数据，精确定位出720张，定位准确率为95.25％，定位一张原始的虹膜图像仅需3-4s(在Matlab下)。能准确定位的图像识别的准确率为97.2％。综上所述，利用本发明的方法可以快速准确的实现相同虹膜和不同虹膜的识别。

步骤1.采集虹膜图像：

这里是利用已经采集好的CASIC 1.0虹膜库进行处理，部分虹膜图像如图8所示：

下面重点分析一个图像029_1_3，见图4(a)，其他的都是使用一样的方法进行处理。

步骤2.显示原始灰度图像I(x，y)的直方图，如图3所示，在左边第一个峰值为瞳孔区域的灰度值，取阈值T0为其向右的第一个谷值，如图3中T0所示：

步骤3.对图4(a)中得到的原始灰度图像I(x，y)，用阈值T0进行二值化处理，得到如图4(b)所示的图像。

步骤4.瞳孔粗定位，由二值化后的图灰度值投影情况，见图4(c)，得到瞳孔粗定位圆心为(178，149)、半径为41。

步骤5.瞳孔精定位得到圆心为(174，149)、半径为42，定位效果如图9所示。

步骤6.虹膜外缘粗定位，得到的Avg＝143.71，双阈值对原灰度图像I(x，y)进行处理并中值滤波后如图4(d)，得到虹膜外缘圆心为(174，149)、半径为109。

步骤7.虹膜外缘精定位，得到虹膜精定位圆心为(177，148)、半径为110.8，虹膜内外圆定位如图4(e)所示。

步骤8.对原图像进行归一化处理，得到如图4(f)所示的图像。

步骤9.对图8得到的图像进行光照不均的校正，得到如图4(g)的图像。

步骤10.对虹膜(图4(g))进行编码，图4(g)的虹膜编码如下：

表2虹膜编码

0010111…

1101101

0000100

1110001

1101111

0000110

0000001

1111110

1101111

0000110

0111001

0010010

0011100

1110100

1000001

0111100

0001000

0101010

1010010

0010010

0000001

1101111

1001111

0010000

1000000

0111110

0000000

1111000

0000000

0111100

1010000

0001000

…………………

可以用图片表示为二值图像，如图6(a)所示，其中白色点为虹膜特征点。

步骤11.虹膜匹配，虹膜029_1_3和029_1_2的虹膜特征图像分别如图6(a)、(b)，

求得它们的Hamming距离为0.246，而由整个虹膜库得到的判定阈值为0.400，由于0.246＜0.400，所以可以判定他们为同一个虹膜；否则如果Hamming距离值大于0.400，可以判定为不同的虹膜。

步骤12.对于整个虹膜库进行处理，对于一个虹膜图像，与之相同的类之间两两比较，不同类之间比较这个类的第一幅虹膜图像，得到了相同类之间和不同类之间的Hamming距离分布图如图7所示。

计算得到FAR为1.5％、FRR为4.3％、EER2.9％、识别准确率为97.1％

Claims (1)

1.一种基于小波变换过零表示的虹膜身份识别算法，其特征在于包括如下步骤：

步骤1，将虹膜摄像采集装置靠近人眼，在人眼进入摄像范围内后进行图像采集，通过数字图像采集CCD的光电转换得到含有虹膜图像的原始灰度图像I(x，y)，图像左上角为I(x，y)的零点I(0，0)，图像向右令横坐标x递增，图像向下令纵坐标y递增，得到人眼图像，虹膜图像的原始灰度图像I(x，y)的取值范围为0～255，x的范围为1～M，y的范围为1～N，原图像大小为M×N；

步骤2，显示原始灰度图像I(x，y)的直方图，标识1左面的第一个峰值为瞳孔区域的灰度值，取阈值T0为其向右的第一个谷值；

步骤3，对步骤1中得到的原始灰度图像I(x，y)，用阈值T0进行二值化处理；

原始灰度图像I(x，y)小于T0的点I(x，y)的值取为0，I(x，y)大于T0的点I(x，y)的值取为1，组成二值矩阵I’(x，y)，得到的二值化图像I’(x，y)大小等于原图像I(x，y)；

步骤4，瞳孔粗定位；

通过步骤3得到原始灰度图像I(x，y)的二值化后的图像I’(x，y)，用式子 $P\left(x\right)=\underset{y}{\mathrm{\Σ}}{I}^{\′}(x,y),$ $P\left(y\right)=\underset{x}{\mathrm{\Σ}}{I}^{\′}(x,y)$ 将二值化后的图像进行水平方向上和垂直方向上的投影，得到二值化后的图像水平方向上和垂直方向上的投影值P(x)，P(y)；并且通过计算 $X_p=\arg \underset{x}{\min }P\left(x\right),$ $Y_p=\arg \underset{y}{\min }P\left(y\right),$ 得到瞳孔粗定位圆的圆心在直角坐标上的结果(x_p，y_p)，如说明书附图3(c)所示；从水平的投影值P(x)谷值点X_p向左搜索，得到点P1，满足在P1以后的点的投影值变化比较平缓(在5以内)，同理，从谷值点X_p向右搜索得到P2，满足在P2以后的点的投影值变化比较平缓(在5以内)；取P1、P2之间部分的范围为瞳孔的粗定位直径，并且可以得到瞳孔粗定位圆的半径r_p＝P2-P1；

步骤5，瞳孔精定位方法：根据步骤4的瞳孔粗定位结果，圆心(x_p，y_p)、半径r_p，通过迭代运算求得瞳孔精确定位结果，圆心(x₀，y₀)、半径r，方法为：

$\max (r,x_0,y_0)|\underset{m}{\mathrm{\Σ}}{I}^{\′}(x_{m,k},y_{m,k})-\underset{m}{\mathrm{\Σ}}{I}^{\′}(x_{m,k-1},y_{m,k-1})|,---\left(1\right)$

其中X_m，k＝r×cos(m×Δθ)+x₀，x₀＝x_p+α，α＝-10～+10，

y_m，k＝r×sin(m×Δθ)+y₀，y₀＝yp+β，β＝-10～+10，

r＝r_p+k×Δr，k＝-10～+10，α、β、k为变量

Δθ是圆周上的一个角度采样间隔，取Δθ＝1°，所以m＝0～360，Δr为在半径方向上的间隔，取Δr＝1；

通过α、β、k的代入运算，求得当式子(1)得到最大值时的x₀、y₀、r的值，则瞳孔的精确定位结果：圆心为(x₀，y₀)、半径为r₀；

步骤6，虹膜外缘粗定位方法：从步骤5得到的瞳孔精定位结果，由虹膜外缘半径范围的先验知识，取一个与瞳孔同心、半径值为比瞳孔半径值大20个像素的同心圆，求出这个圆和瞳孔圆之间的虹膜部分的灰度值平均值Avg：

Avg＝位于这环形区域点的灰度值的和÷位于这环形区域之间的总点数；

计算过程是：

环形区域中的点的集合s：{(x，y)|瞳孔圆心²≤x²+y²≤(瞳孔圆心+20)²}原图像上的灰度值为I(x，y)，位于这环形区域点的灰度值的和为

位

于这环形区域之间的总点数为集合s的点的总数；

令阈值T1＝Avg-15，阈值T2＝Avg+15，作为双阈值对原灰度图像I(x，y)进行如下处理：

这样得到图像矩阵I”(x，y)是基本上只包含虹膜部分的图像，再对其进行中值滤波，得到图像矩阵I₂(x，y)；

求I₂(x，y)在水平投影值 $P^{\′}\left(x\right)=\underset{v}{\mathrm{\Σ}}{I}_2(x,y),$ 然后求 ${X_p}^{\′}=\arg \underset{x}{\min }{P}^{\′}\left(x\right),$ 从谷值点X_p′向左搜索，得到点P1’，满足在P1’以后的点的点间灰度值变化比较平缓(在5以内)；同理，从谷值点X_p′向右搜索得到P2’，满足在P2’以后的点的点间灰度值变化比较平缓(在5以内)；取P1’、P2’之间部分的范围为虹膜外缘的粗定位直径，并且可以得到虹膜外缘粗定位圆的半径R’＝P2’-P1’，虹膜外缘粗定位圆心设为与瞳孔精定位圆心一样，为(x₀，y₀；

步骤7，虹膜外缘精定位方法：对步骤3二值化后的图像，进行基于Canny算子的边缘提取，得到关于边缘的二值图像，1为边缘，0为其他点；根据步骤6的虹膜外缘粗定位半径R’和圆心(x₀，y₀)，在R’±δ，x₀±α，y₀±β范围内和在图像的-45°～45°and 135°～225°范围内进行Hough变换，其中δ、α、β为待定常数，其取值与精度和运算量有关，它们的值越大精度越高、运算量越大，反之亦然；这样就可以得到虹膜外缘的精定位结果：半径为R，圆心为(x_i，y_i)；

步骤8，对原图像进行归一化处理：由步骤5的瞳孔精定位结果(x₀，y₀，r₀)和步骤7的虹膜外缘精定位结果(x_i，y_i，R)，进行虹膜图像归一化处理；把原灰度图像的虹膜部分归一化为64×512的矩形图像I_nor(i，j)，1≤i≤64、1≤j≤512；

步骤9，对步骤8得到的图像进行光照不均的校正，然后进行直方图均衡化，得到总体亮度均匀、纹理清晰的图像I_adj(i，j)，1≤i≤64、1≤j≤512；

步骤10，虹膜编码：对步骤9所得到的图像进行特征点提取；先用4阶B-样条小波滤波器对步骤9所得到的图像I_adj(i，j)分别在水平方向和垂直方向进行卷积滤波；4阶B-样条小波滤波器：是一个小波滤波器，其参数如下：

表1 4阶B-样条小波滤波器

K是变量，ψ(K)是相应的小波函数值；卷积方法如下：

temp_h(i，j)＝ψ(-i)*I_adj(i，j)……行方向上卷积，j固定

temp_v(i，j)＝ψ(-j)*I_adj(i，j)……列方向上卷积，i固定

分别得到行和列方向上的卷积滤波后的临时存储矩阵temp_h(i，j)，temp_v(i，j)，并且进行如下判决得到行和列方向上编码，行和列方向上的编码进行与操作得到虹膜编码：

3.令code(i，j)＝code_h(i，j)And code_v(i，j)，这里的And为求两个数的与运算，以上i＝1～64，j＝1～512；

这样二值化矩阵code(i，j)就是由虹膜特征点所组成的虹膜编码，大小为64×512；

步骤11，编码匹配：由步骤9得到的矩阵code(i，j)可以代表原始图像的虹膜的编码，如果要把两个虹膜编码code1(i，j)和code2(i，j)进行匹配，可以使用下面的Hamming距离的式子，其表征的意义为两个虹膜编码之间的相似程度：

$\mathrm{Ham}\min g\_\mathrm{Dis}\tan \mathrm{ce}=\underset{\mathrm{\φ}}{\mathrm{Min}}\left\{\frac{1}{M\×N}\underset{i=0}{\overset{M}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}\left(\mathrm{code}1(i+\mathrm{\φ},j)\right)\mathrm{Xor}\left(\mathrm{code}2(i,j)\right)\right\}$

for-10≤φ≤10

如果HammingDistance小于某阈值T，则两个虹膜是相同虹膜；如果HammingDistance大于T，则两个虹膜是不相同虹膜；其中T的取值和样本的性质有关，可以在数据统计的时候决定最优的阈值T。

Images