CN100351851C - 基于数学形态学和概率统计的虹膜定位方法 - Google Patents

Abstract

本发明提供一种基于数学形态学和概率统计相结合的虹膜定位方法。它首先二值化虹膜图像，利用投影得到虹膜内圆粗定位圆心，然后用运用形态学进行边缘提取得到虹膜内圆的边缘，以内圆粗定位圆心(及其周围的点)到其边缘上每点的距离作为一个样本，从而将搜索圆心和半径三个参数转换为搜索最小均值的样本，而此样本的均值就是圆的半径。定出内圆后，用二值化和数学形态学方法得到外圆的边缘，然后以内圆圆心为外圆粗定位圆心，用定位内圆的方法定位外圆。本发明提供虹膜定位方法与现有的虹膜定位方法相比，降低了标准误差，提高了精确度和速度。

Description

基于数学形态学和概率统计的虹膜定位方法

技术领域

基于数学形态学和概率统计的虹膜定位方法，属于图像信号处理技术领域，特别涉及虹膜识别技术中的虹膜定位方法。

背景技术

生物识别是当前信息安全发展的主要方向，是目前世界信息安全领域的前沿研究课题。虹膜生物识别技术则是目前生物识别技术中易操作、精度高、最有市场前景的技术。它是解决信息安全问题的重要手段，是计算机图像处理技术和模式识别技术在信息安全领域的应用。该技术在解决信息安全领域的诸多问题，如门禁控制系统，网络技术的信息安全，网络经济的电子商务，个人身份的精确认证等具体问题中应用广泛，具有巨大的经济价值和现实意义。同时，它在身份证技术的发展和社会安全管理中占有重要地位。随着计算机性能的提高和图像处理技术的发展，虹膜生物识别系统已经日趋成熟。详见文献：Clarke R.HumanIdentification System：Management Challenges and Public Policy Issues.InformationTechnology&People，1994，7(4)：6～37和文献：范君.虹膜识别生物技术激发10亿安全市场.金卡工程，2003，5：61-64所述。

在虹膜生物识别技术中，虹膜的定位是整个识别技术的重点和难点。虹膜的定位时间和精度直接影响整个虹膜生物识别系统的性能。虹膜定位是指对人的眼睛的虹膜部分进行精确地定位，这是虹膜识别至关重要而关键的一步，虹膜能否准确地定位，关系到下一步虹膜纹理特征的能否准确地提取，从而对分析结果产生重要的影响。而在虹膜定位中，虹膜外边缘的提取则是一个难点。因此，如何在现有所有有价值的研究成果之上，缩短定位时间和提高定位精度特别是外边缘的定位精度将成为我们当前研究的主要方向。详见文献：王春，叶虎年.虹膜识别算法的研究[J].贵州工业大学学报(自然科学版)，2000，29(3)：48 52.和文献：李庆嵘，马争.虹膜定位算法研究[J].电子科技大学学报，2002，31(1)：79.和文献：Daugman.Statisticalrichness of visual phase information Update on recognizing persons by their iris pattern.International Journal ofcomputer Vision，2001，45(1)：25-38.所述。

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

(1)基于边缘提取和Hough变换的虹膜定位算法。它通过边缘检测算子来提取虹膜的内外边缘，然后运用Hough变换来定位虹膜的内外圆。其缺点是Hough变换需要设置阈值。由于受光线的影响，瞳孔的大小是变化的，因而虹膜边缘上点的数量也是变化的，对于大量的虹膜图像而言，一个准确而合适的阈值是很难得到的。而且，该方法随点数的增加，运算量增大速度变慢。详见文献：王春，叶虎年.虹膜识别算法的研究[J].贵州工业大学学报(自然科学版)，2000，29(3)：4852.

(2)两步法虹膜定位算法。它首先根据眼部的灰度分布特点做x方向灰度投影和y方向灰度投影进行粗定位，然后运用算子

进行精确定位。其缺点积分运算量大，速度慢。详见文献：Daugman.Statistical richness of visual phase informationUpdate on recognizing persons by their iris pattern.International Journal of computer Vision，2001，45(1)：25-38.

发明内容

本发明的任务是提供一种基于数学形态学和概率统计的虹膜定位方法，它具有定位精度高和速度快等特点。

为了方便地描述本发明内容，首先介绍几个概念，并对一些术语进行定义。

概念一：数学形态学。数学形态学以图像分析为基础，用具有一定形态结构的“结构元素”去度量图像的形态，以解决图像理解问题。形态学的基础是腐蚀和膨胀运算，以及由此而产生的开、闭运算。腐蚀和膨胀运算的公式分别为：U＝AΘB＝{U：B+UA}，V＝AB＝{V：(-B+V)∩A≠Ф}；开和闭运算的公式分别为：AοB＝(AΘB)B，A·B＝(AB)ΘB。其中，A为原始图像，B为结构元素，U为原始图像经腐蚀后得到的图像，V为原始图像经膨胀后得到的图像，Θ为腐蚀运算符号，为膨胀运算符号，о为开运算符号，·为闭运算符号，Ф为空集符号。

概念二：概率统计。假设数组{x}服从均值是u＝2.5，方差是σ²的正态分布。根据概率统计的理论，随着σ²的减小，数组{x_i}与均值的差值的绝对值和也就相应的减小。也就是说，数组{x}的值的分布随着σ²的减小而越来越集中于均值的两侧。因此，当σ²无限趋近于0时，数组{x_i}的值的分布也就无限趋近于均值。而若一个圆上所有点到圆心的距离作为一数组{x_i}，则数组{x_i}的值均为半径r，标准差为0。本发明以圆内某点到其边缘上点的距离作为一个样本值。以圆内某点到其边缘上所有点的距离作为一个样本。

概念三：标准误差。设n个测量值的误差为ε₁、ε₂......ε_n，则这组测量值的标准误差σ等于： $\mathrm{\σ}=\sqrt{({\mathrm{\ϵ}}_1^2+{\mathrm{\ϵ}}_2^2+\mathrm{\Λ\Λ}+{\mathrm{\ϵ}}_n^2)/n}=\sqrt{\mathrm{\Σ}{\mathrm{\ϵ}}_1^2/n}.$ 由于被测量的真值是未知数，各测量值的误差也都不知道，因此不能按上式求得标准误差。测量时能够得到的是算术平均值( N)，它最接近真值(N)，而且也容易算出测量值和算术平均值之差，称为残差(记为v)。理论分析表明可以用残差v表示有限次(n次)观测中的某一次测量结果的标准误差σ，其计算公式为： $\mathrm{\σ}=\sqrt{[{(N_1-N)}^2+{(N_2-N)}^2+\mathrm{\Λ\Λ}{(N_n-N)}^2]/n}=\sqrt{\mathrm{\Σ}{\mathrm{\υ}}_i^2/n}$ 需要注意的是，标准误差不是测量值的实际误差，也不是误差范围，它只是对一组测量数据可靠性的估计。标准误差小，测量的可靠性大一些，反之，测量就不大可靠。

定义1：虹膜。位于人眼中巩膜和瞳孔之间的部分称为虹膜。不同人之间的虹膜具有随机的细节特征和纹理图像，这些特征在人的一生中保持相当的稳定性，它受眼皮、角膜的保护不易改变。

定义2：二维中值滤波：二维中值滤波(Median filtering)是基于排序统计理论的一种能有效抑制噪声的非线性信号处理技术，它在滤除噪声(尤其是脉冲噪声)的同时能很好地保护信号的细节信息。其基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代换。

定义3：二值化过程。把整幅图像的所有值转化成只有两种值的过程，一般这两种值为0和1或者0和255。其具体的方法是：当图像上的值大于或者等于二值化阀值的时候，该点的值被二值化为1(或255)；当图像上的值小于二值化阀值的时候，该点的值被二值化为0。

定义4：灰度直方图。图像的灰度级范围为0，1，...，255，设灰度级i的像素数为n_i，则一幅图像的总像素为 $N=\underset{i=0}{\overset{255}{\mathrm{\Σ}}}{n}_i,$ 灰度级i出现的概率定义为 $P_i=\frac{n_i}{N}.$ 灰度直方图即为灰度级i的像素数n_i与灰度级i的二维关系，它反映了一幅图像上灰度分布的统计特性，成为利用像素灰度作属性的分割方法的基础。

定义5：水平投影。一种通过水平方向累加，现实把二维空间中的图像灰度水平分布值转换到一维空间中的方法，该变换函数为 $P\left(x\right)=\underset{y}{\mathrm{\Σ}}I(x,y).$ 其中I(x，y)表示第x行第y列的图像灰度值，P(x)为第x行的水平投影值，y的变化范围是从1到n，n代表图像的长度。

定义6：垂直投影。一种通过垂直方向累加，现实把二维空间中的图像灰度垂直分布值转换到一维空间中的方法，该变换函数为 $P\left(y\right)=\underset{x}{\mathrm{\Σ}}I(x,y).$ 其中I(x，y)表示第x行第y列图像灰度值，P(y)为第y列的垂直投影值，x的变化范围是从1到m，m代表图像的宽度。

定义7：形态学边缘提取。集合A的边界记为β(A)，可以通过下述算法提取边缘：设B是一个合适的结构元素，首先令A被B腐蚀，然后求集合A和它的腐蚀的差。如下式所示：β(A)＝A-(AΘB)。

定义8：外圆二值化阈值。对虹膜图像外圆进行二值化时所选用的门限。其计算公式为： $\mathrm{yuzhi}=\left[\frac{1}{2}(\mathrm{aver}1+\mathrm{aver}2)\right]/256,$ 其中yuzhi为二值化阀值，aver1为经过二维中值滤波后的圆环的均值(圆环以内圆圆心为圆心，圆环内圆半径为内圆半径+5，圆环外圆半径为内圆半径+10)，aver2是两个11×11大小的方块的均值(方块的中心分别是以内圆圆心的行坐标为行坐标，以内圆圆心的列坐标减130为列坐标和以内圆圆心的行坐标为行坐标，以内圆圆心的列坐标加140为列坐标)，/为除号，256为归一化yuzhi。

定义9：维纳滤波。维纳滤波器是一种线性平滑滤波器，它是一种自适应滤波器，能根据图像的区域方差来调整滤波器的输出 ${\mathrm{\δ}}^2=\frac{1}{L}\underset{(x,y)\∈l}{\mathrm{\Σ}}[F^2(x,y)-G^2(x,y)]$ 其中，L为选取的邻域M×N，F(x，y)调整前的像素点值，G(x，y)是邻域平均值(如定义5所示)，δ为均方差。通过邻域m×n估算平均值和标准偏差，对图像应用像素平滑自适应滤波： $F^{*}(x,y)=G(x,y)+\frac{{\mathrm{\δ}}^2-v^2}{{\mathrm{\δ}}^2}[F(x,y)-G(x,y)],$ 其中，v²为噪声方差。

本发明详细技术方案为：

基于数学形态学和概率统计的虹膜定位方法，它包含下列步骤：

步骤1、通过虹膜图像摄像装置(为了获得能够进行纹理分析的虹膜图像而专门制造的一种图像摄像装置)，获得能够进行虹膜纹理分析的眼部图像。该图像是只包含亮度信息而没有任何其他颜色信息的灰度图像。

步骤2、对输入的虹膜图像进行7×7窗口的二维中值滤波，达到平滑虹膜图像的目的，减弱光照的不均匀对二值化的影响。

步骤3、虹膜内圆圆心的粗定位，具体包括以下步骤：

步骤1)、求出二维中值滤波后的虹膜图像的灰度直方图的第一个波峰的极值，以该取得极值时的灰度级i+5为门限，如果I(x，y)＞i+5，令I(x，y)＝0；如果I(x，y)＜i+5，令I(x，y)＝1，其中I(x，y)表示第x行第y列的图像强度值对虹膜图像进行二值化。然后用结构元素se(se是半径为15的圆盘)对二值图像进行形态开运算，剔除瞳孔部分的其他孤立点和面，平滑虹膜内圆边缘。

步骤2)、对经过形态开运算处理的二值图像进行水平和垂直投影得到虹膜内圆的粗定位运用式 $P\left(x\right)=\underset{y}{\mathrm{\Σ}}I(x,y)$ (其中I(x，y)表示第x行第y列的图像强度值，P(x)为第x行的水平投影值，y的变化范围是从1到n，n代表图像的长度)求出二值图像每行的水平投影值，以水平投影值最大的那行的行坐标作为粗定位内圆圆心的行坐标a₁₁；运用式 $P\left(y\right)=\underset{x}{\mathrm{\Σ}}I(x,y)$ (其中I(x，y)表示第x行第y列的图像强度值，P(y)为第y列的垂直投影值，x的变化范围是从1到m，m代表图像的宽度)求出二值图像每列的垂直投影值，以垂直投影值最大的那列的列坐标作为粗定位内圆圆心的列坐标b₁₁。

步骤4、虹膜内圆的精确定位，具体包括以下步骤：

步骤1)、对经过形态开运算处理的二值图像进行形态学边缘提取

运用式β(A)＝A-(AΘB)(其中β(A)为集合A的边界，B是一个合适的结构元素)对经过步骤3处理的二值图像进行边缘提取，得到虹膜的内边缘。

步骤2)、以虹膜内圆的粗定位圆心(a₁₁，b₁₁)为基点，计算虹膜的内边缘上每点到点(a₁₁，b₁₁)的距离 ${\stackrel{\mathrm{\ω}}{D}}_1^2={(\stackrel{\mathrm{\ω}}{i}-a_{11}\·E)}^2+{(\stackrel{\mathrm{\ω}}{j}-b_{11}\·E)}^2$ (其中

表示内边缘上每点的坐标，E为与向量 同大小的值均为1的向量，

表示一个样本)。用式 ${\stackrel{\mathrm{\ω}}{D}}_2^2={(\stackrel{\mathrm{\ω}}{i}-a_{11}\·E)}^2+{(\stackrel{\mathrm{\ω}}{j}-(b_{11}+1)\·E)}^2$ (其中

表示内边缘上每点的坐标，E为与向量

同大小的值均为1的向量，

表示一个样本)计算虹膜内边缘上每点到点(a₁₁，b₁₁+1)的距离。继续计算虹膜内边缘上每点到点(a₁₁，b₁₁+2)的距离

直到计算出虹膜内边缘上每点到点(a₁₁+9，b₁₁+9)的距离

步骤3)、根据公式 ${\mathrm{\σ}}^2=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-u)}^2$ (其中x_i表示样本的样本值，n表示样本的容量，u表示样本均值)，计算每个样本

(k＝[1,100])的方差，找出方差最小的样本 根据公式 $u=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$ (其中x_i表示样本的样本值，n表示样本的容量)，求出样本

的均值r；根据z的值求出对应的点的坐标

(其中qs表示求整数商，qy表示求整数余数)。则点的坐标

即为内圆的圆心(a₁，b₁)，均值r为内圆半径。

步骤5、虹膜外圆的定位，具体包括以下步骤：

步骤1)、根据公式 $\mathrm{yuzhi}=\left[\frac{1}{2}(\mathrm{aver}1+\mathrm{aver}2)\right]/256$ (其中yuzhi为二值化阀值，aver1为经过二维中值滤波后的圆环的均值(圆环以内圆圆心为圆心，圆环内圆半径为内圆半径+5，圆环外圆半径为内圆半径+10)，aver2是两个11×11大小的方块的均值(方块的中心分别是以内圆圆心的行坐标为行坐标、以内圆圆心的列坐标减130为列坐标和以内圆圆心的行坐标为行坐标、以内圆圆心的列坐标加140为列坐标)，/为除号)，计算外圆二值化的阈值。然后以该阈值对二维中值滤波后的虹膜图像进行二值化，得到图像BW1。

步骤2)、用结构元素se(se是半径为15的圆盘)对图像BW1进行形态开运算，得到图像BW2，从而得到图像BW3＝BW1-BW2。对图像BW1进行补运算得到图像BW4，然后得到图像BW5＝BW4+BW3。对图像BW5进行区域填充，得到图像BW6。填充二值图像中的孔洞，消除虹膜外边缘内的孔洞对于虹膜外边缘提取的影响。对图像BW6进行形态学边缘提取。运用式β(A)＝A-(AΘB)(其中β(A)为集合A的边界，B是一个合适的结构元素)对图像BW6进行边缘提取，得到虹膜的外边缘。由于受上眼皮和下眼皮的影响，因此仅仅取左右两边的一部分圆弧(根据内圆圆心和虹膜的外边缘来获取)来进行定位。

步骤3)、以虹膜内圆圆心(a₁，b₁)为基点，计算虹膜的外边缘上每点到点(a₁，b₁)的距离 ${\stackrel{\→}{D}}_1^{\′2}={(\stackrel{\→}{i}-a_1\·E)}^2+{(\stackrel{\→}{j}-b_1\·E)}^2$ (其中 表示外边缘上每点的坐标，E为与向量 同大小的值均为1的向量，

表示一个样本)。用式 ${\stackrel{\→}{D}}_2^{\′2}={(\stackrel{\→}{i}-a_1\·E)}^2+{(\stackrel{\→}{j}-(b_1+1)\·E)}^2$ (其中

表示外边缘上每点的坐标，E为与向量

同大小的值均为1的向量， 表示一个样本)，计算虹膜外边缘上每点到点(a₁，b₁+1)的距离。继续计算虹膜内边缘上每点到点(a₁，b₁+2)的距离

直到计算出虹膜内边缘上每点到点(a₁+9，b₁+9)的距离

步骤4)、对每个样本进行维纳滤波，以减少非虹膜外边缘点对定位的影响。然后根据公式 ${\mathrm{\σ}}^2=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-u)}^2$ (其中x_i表示样本的样本值，n表示样本的容量，u表示样本均值)，计算每个样本 ${\stackrel{\→}{D}}_k^{\′}(k=[\mathrm{1,100}\left]\right)$ 的方差，找出方差最小的样本

根据公式 $u=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$ (其中x_i表示样本的样本值，n表示样本的容量)，求出样本

的均值R；根据z的值求出对应的点的坐标

(其中qs表示求整数商，qy表示求整数余数)。则点的坐标

即为外圆的圆心(a₂，b₂)，均值R为外圆半径。

通过以上步骤，就可以从含有虹膜的原始图像中提取出虹膜。

需要说明的是：

1.步骤1中的虹膜图像来源于中科院自动化所的CASIA1.0虹膜数据库。

2.步骤3的分步骤(1)中用结构元素se(se为圆盘，半径为15)进行形态开运算，而不用其他结构元素，是因为实验结果表明二值化图像其虹膜内圆外的连通对象均小于半径为15的圆盘，用半径为15的圆盘结构元素能够完全消除内圆外的其他点的影响。

3.在步骤3的分步骤(2)中，实际上水平和垂直投影的最大值都不止一个，而是相邻的好几个。这是因为在计算机中图像数据是以矩阵的形式表示的，而用矩阵表示圆则在圆周上会出现很短的线段。而步骤3的分步骤(1)处理后的内圆并不是一个规则的圆盘，仿真表明其水平和垂直投影的最大值是相邻的一段。本发明选取最先出现的最大值的坐标加上一个偏移量，作为粗定位的内圆圆心，而在实际中，根据不同的虹膜库来调整偏移量。

4.步骤4中的分步骤(1)和步骤5中的分步骤(2)中，运用形态学进行边缘提取是因为运用形态学进行边缘提取比运用边缘检测算子进行边缘提取的速度要快。

5.步骤5中的分步骤(2)中，对步骤5中的分步骤(1)中的二值图像进行开运算、图像间的减法运算和图像间的加法运算，是为了进一步的提取虹膜外边缘。而进行区域填充的目的则是消除为虹膜外边缘内的孔洞，使外边缘内的区域完全连通，以使得形态学边缘提取时提取的边缘仅仅时虹膜的外边缘，而没有其他点的影响。而仅仅取左右两边的一部分圆弧(根据内圆圆心和估计的外圆半径来获取)来进行定位，是由于大多数虹膜图像其上眼皮或下眼皮经常是遮盖了虹膜外边缘的上边界或下边界，而进行虹膜外圆定位仅需要左右两边的圆弧就可以精确定位。因此，为了算法的通用性和避免不必要的计算，仅仅取左右两边的一部分圆弧来进行定位。

6.步骤5的分步骤(4)中，对每个样本进行维纳滤波，是为了平滑样本值，减小非边缘点的样本值对外边缘定位的影响。

7.步骤4的分步骤(2)和步骤5的分步骤(3)中的搜索范围是根据CASIA1.0虹膜数据库中的图像仿真实验后得到的。而对于不同的虹膜库可以适当的扩大搜索的范围。

本发明提出的基于数学形态学和概率统计相结合的方法，可以有效地提高定位的速度和定位精度等性能。运用形态开运算能够平滑图像的轮廓，削弱狭窄的部分，去掉细的突出。本发明中形态开运算可以消除虹膜内圆之外和虹膜外圆之外的小于结构元素的非边界部分对虹膜定位的影响。因而对步骤3的二值图像进行形态开运算可以使得其仅仅剩余瞳孔部分，对步骤13的二值图像进行形态开运算和图像的加、减法运算可以使得其仅仅剩余瞳孔部分、虹膜部分和上眼皮部分或下眼皮部分。

本发明的方法采用数学形态学和概率统计相结合的方法，首先采用中值滤波、二值化和数学形态学的方法突出了虹膜的内圆区域和外圆区域，实现了虹膜的内外边缘的提取；然后运用标准误差原理将每次定位得到的圆看作为有限次(n次)观测中的某一次测量结果，而测量结果的标准误差σ则是对一组测量数据可靠性的估计，标准误差小，测量的可靠性大一些，反之，测量就不大可靠，再根据概率统计的方法实现了将虹膜圆心和半径的定位转化为样本的方差和均值的求取。采用本发明的虹膜定位方法，不仅可以准确定位虹膜的内外边缘，而且方法的计算速度也满足了实时性的需要。

本发明的创新之处在于：

1、结合形态学相关知识，利用中值滤波和通过合理选择二值化阀值对中值滤波后的图像进行二值化的方法，在原始图像中突出虹膜的内圆和外圆位置，实现了虹膜内外边缘的提取，特别是外边缘的提取。本发明的提取外边缘的方法不仅迅速，而且很精确。

2、运用标准误差原理，将每次定位得到的圆看作为有限次(n次)观测中的某一次测量结果，而测量结果的标准误差σ则是对一组测量数据可靠性的估计，标准误差小，测量的可靠性大一些，反之，测量就不大可靠。根据概率统计的方法实现了将虹膜圆心和半径的定位转化为样本的方差和均值的求取。

3、基于数学形态学和概率统计相结合的方法，提高了虹膜定位的速度，满足了虹膜识别系统实时性的需要。首先利用二值化和形态学的方法，提取出虹膜的内外边缘，然后运用标准误差原理和概率统计原理再对虹膜的内外边缘进行精确定位。

附图说明

图1是含有虹膜的原始虹膜图像示意图；

其中，1表示人的上眼皮；2所在的圆环部分表示人的虹膜；3表示人的瞳孔；4表示人的巩膜；5表示人的下眼皮。

图2是本发明最终得到的虹膜定位图像示意图；

其中，大圆代表虹膜的外边缘，小圆代表虹膜的内边缘。

图3是本发明的流程示意图。

具体实施方法

采用本发明的方法，首先使用Matlab语言编写虹膜定位软件，然后中科院自动化所的CASIA1.0虹膜数据库为源数据输入到膜定位软件中进行处理；经过虹膜内圆定位和外圆定位后得到一幅包含可进行纹理分析的虹膜图像。采用共812张不同人的和相同人的不同时期的虹膜图像作为源数据，准确定位出800张，定位的准确率为98.52％，定位一幅虹膜图像所需时间不超过1s。

综上所述，本发明的方法充分利用形态学和概率统计的特性，结合虹膜眼部图像的灰度分布特征，从而实现快速准确地从所提供的虹膜图像中定位出虹膜。

Claims (1)

1、基于数学形态学和概率统计的虹膜定位方法，其特征在于，它包含下列步骤：

步骤1、通过虹膜图像摄像装置获得能够进行虹膜纹理分析的眼部图像；

步骤2、对输入的虹膜图像进行7×7窗口的二维中值滤波，达到平滑虹膜图像的目的，减弱光照的不均匀对二值化的影响；

步骤3、虹膜内圆圆心的粗定位，具体包括以下步骤：

步骤1)、求出二维中值滤波后的虹膜图像的灰度直方图的第一个波峰的极值，以该取得极值时的灰度级i+5为门限，如果I(x，y)＞i+5，令I(x，y)＝0；如果I(x，y)＜i+5，令I(x，y)＝1，其中I(x，y)表示第x行第y列的图像强度值，对虹膜图像进行二值化；然后用半径为15的圆盘结构元素se对二值图像进行形态开运算，剔除瞳孔部分的其他孤立点和面，平滑虹膜内圆边缘；

步骤2)、对经过形态开运算处理的二值图像进行水平和垂直投影得到虹膜内圆的粗定位

运用式 $P\left(x\right)=\underset{y}{\mathrm{\Σ}}I(x,y),$ 其中I(x，y)表示第x行第y列的图像强度值，P(x)为第x行的水平投影值，y的变化范围是从1到n，n代表图像的长度，求出二值图像每行的水平投影值，以水平投影值最大的那行的行坐标作为粗定位内圆圆心的行坐标a₁₁；运用式 $P\left(y\right)=\underset{x}{\mathrm{\Σ}}I(x,y)$ 其中I(x，y)表示第x行第y列的图像强度值，P(y)为第y列的垂直投影值，x的变化范围是从1到m，m代表图像的宽度，求出二值图像每列的垂直投影值，以垂直投影值最大的那列的列坐标作为粗定位内圆圆心的列坐标b₁₁；

步骤4、虹膜内圆的精确定位，具体包括以下步骤：

步骤1)、对经过形态开运算处理的二值图像进行形态学边缘提取

运用式β(A)＝A-(AΘB)，其中β(A)为集合A的边界，B是个合适的结构元素，对经过步骤3处理的二值图像进行边缘提取，得到虹膜的内边缘；

步骤2)、以虹膜内圆的粗定位圆心(a₁₁，b₁₁)为基点，计算虹膜的内边缘上每点到点(a₁₁，b₁₁)的距离 ${\stackrel{\→}{D}}_1^2={(\stackrel{\→}{i}-a_{11}\·E)}^2+{(\stackrel{\→}{j}-b_{11}\·E)}^2,$ 其中

表示内边缘上每点的坐标，E为与向量

同大小的值均为1的向量，

表示一个样本；用式 ${\stackrel{\→}{D}}_2^2={(\stackrel{\→}{i}-a_{11}\·E)}^2+{(\stackrel{\→}{j}-(b_{11}+1)\·E)}^2,$ 其中

表示内边缘上每点的坐标，E为与向量

同大小的值均为1的向量，

表示一个样本，计算虹膜内边缘上每点到点(a₁₁，b₁₁+1)的距离；继续计算虹膜内边缘上每点到点(a₁₁，b₁₁+2)的距离直到计算出虹膜内边缘上每点到点(a₁₁+9，b₁₁+9)的距离

步骤3)、根据公式 ${\mathrm{\σ}}^2=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-u)}^2,$ 其中x_i表示样本的样本值，n表示样本的容量，u表示样本均值，计算每个样本

(k＝[1，100])的方差，找出方差最小的样本

根据公式 $u=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i,$ 其中x_i表示样本的样本值，n表示样本的容量，求出样本 的均值r；根据z的值求出对应的点的坐标 $(a_{11}+\mathrm{qs}\left(\frac{z}{10}\right),b_{11}+\mathrm{qy}\left(\frac{z}{10}\right)),$ 其中qs表示求整数商，qy表示求整数余数；则点的坐标 $(a_{11}+\mathrm{qs}\left(\frac{z}{10}\right),b_{11}+\mathrm{qy}\left(\frac{z}{10}\right))$ 即为内圆的圆心(a₁，b₁)，均值r为内圆半径；

步骤5、虹膜外圆的定位，具体包括以下步骤：

步骤1)、根据公式 $\mathrm{yuzhi}=\left[\frac{1}{2}(\mathrm{aver}1+\mathrm{aver}2)\right]/256,$ 其中yuzhi为二值化阀值，aver1为经过二维中值滤波后的圆环的均值，圆环以内圆圆心为圆心，圆环内圆半径为内圆半径+5，圆环外圆半径为内圆半径+10，aver2是两个11×11大小的方块的均值，方块的中心分别是以内圆圆心的行坐标为行坐标、以内圆圆心的列坐标减130为列坐标和以内圆圆心的行坐标为行坐标、以内圆圆心的列坐标加140为列坐标，/为除号，计算外圆二值化的阈值；然后以该阈值对二维中值滤波后的虹膜图像进行二值化，得到图像BW1；

步骤2)、用半径为15的圆盘结构元素se对图像BW1进行形态开运算，得到图像BW2，从而得到图像BW3＝BW1-BW2；对图像BW1进行补运算得到图像BW4，然后得到图像BW5＝BW4+BW3；对图像BW5进行区域填充，得到图像BW6；运用式β(A)＝A-(AΘB)，其中β(A)为集合A的边界，B是一个合适的结构元素，对图像BW6进行边缘提取，得到虹膜的外边缘；

步骤3)、以虹膜内圆圆心(a₁，b₁)为基点，计算虹膜的外边缘上每点到点(a₁，b₁)的距离 ${\stackrel{\→}{D}}_1^{\′2}={(\stackrel{\→}{i}-a_1\·E)}^2+{(\stackrel{\→}{j}-b_1\·E)}^2,$ 其中

表示外边缘上每点的坐标，E为与向量

同大小的值均为1的向量，

表示一个样本；用式 ${\stackrel{\→}{D}}_2^{\′2}={(\stackrel{\→}{i}-a_1\·E)}^2+{(\stackrel{\→}{j}-(b_1+1)\·E)}^2,$ 其中

表示外边缘上每点的坐标，E为与向量

同大小的值均为1的向量，

表示一个样本，计算虹膜外边缘上每点到点(a₁，b₁+1)的距离；继续计算虹膜内边缘上每点到点(a₁，b₁+2)的距离

直到计算出虹膜内边缘上每点到点(a₁+9，b₁+9)的距离

步骤4)、对每个样本进行维纳滤波，以减少非虹膜外边缘点对定位的影响；然后根据公式 ${\mathrm{\σ}}^2=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-u)}^2,$ 其中x_i表示样本的样本值，n表示样本的容量，u表示样本均值，计算每个样本

(k＝[1，100])的方差，找出方差最小的样本 根据公式 $u=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i,$ 其中x_i表示样本的样本值，n表示样本的容量，求出样本

的均值R；根据z的值求出对应的点的坐标 $(a_1+\mathrm{qs}\left(\frac{z}{10}\right),b_1+\mathrm{qy}\left(\frac{z}{10}\right)),$ 其中qs表示求整数商，qy表示求整数余数；则点的坐标 $(a_1+\mathrm{qs}\left(\frac{z}{10}\right),b_1+\mathrm{qy}\left(\frac{z}{10}\right))$ 即为外圆的圆心(a₂，b₂)，均值R为外圆半径。

Images