CN103198484B - 基于非线性尺度空间的虹膜图像分割算法 - Google Patents

Abstract

本发明具体公开了一种基于非线性尺度空间的虹膜图像分割算法，包括以下步骤：步骤1对虹膜内边界定位：利用非线性尺度空间所具有的尺度演化特性，对原始虹膜图像进行平滑处理，消除瞳孔区域的亮斑，保证瞳孔区域的灰度处于整幅图像的最低水平，再通过阈值实现瞳孔检测，实现虹膜内边界定位；步骤2对虹膜外边界定位：利用非线性尺度空间的尺度演化特性，对原始虹膜图像进行平滑处理，消除几何尺度较小的睫毛对虹膜区域的遮挡，实现了对外边界的准确定位。本算法在外边界定位过程中，同样利用非线性尺度空间的尺度演化特，对原始虹膜图像进行平滑，消除了几何尺度较小的睫毛对虹膜区域的遮挡，实现了对外边界的准确定位。

Description

基于非线性尺度空间的虹膜图像分割算法

技术领域

本发明涉及一种基于非线性尺度空间的虹膜图像分割算法。

背景技术

生物特征识别算法(Biometrics)通过计算机利用人类自身的生理或行为特征进行个人身份认定。这些身体特征包括指纹、虹膜、掌纹、手形、人脸、声音、视网膜和DNA等人体的生理特征，以及签名动作、行走步态、敲击键盘的力度等行为特征。这些生理或行为特征由于人各有异、几乎随身携带并且具有相当的稳定性，能够成为身份鉴定的可靠依据。目前，生物特征识别算法凭借这些独特优势，已经在信息安全、金融交易、社会安全、人员管理、医疗卫生等领域获得了广泛应用。

其他生物特相比，虹膜唯一性高、精度高、稳定性好，适用人群广泛。作为一项重要的生物特征识别算法，虹膜识别得到了科学界和工业界的高度关注。目前，一个典型的虹膜识别系统主要由虹膜图像预处理、特征提取和编码组成。其中虹膜图像预处理的主要任务是正确分割图像中的虹膜区域，即虹膜边界的定位，从而为后续的特征提取和编码过程提供有效信息。一般情况下，人们将虹膜的形状建模为圆环，其内外边界为同心圆，内圆即为瞳孔与虹膜的边界。虹膜预处理的基本步骤为：(1)确定瞳孔中心；(2)确定虹膜内边界；(3)确定虹膜外边界。虹膜边界定位是虹膜识别系统的基础，它会对最终的识别效果产生重要影响。最早的虹膜边界定位方法由Daugman提出，即经典的积分/微分算子(Integro-differentialoperator)；Wildes提出的算法将边缘检测与Hough变换相结合；Boles等采用一维三次样条小波提取出的图像过零点作为特征。随后，人们在虹膜图像预处理阶段引入了多种工具，如活动轮廓模型、傅里叶变换与几何投影、统计方法、Gabor滤波器等。

近年来，考虑到虹膜识别系统的实用性，人们开始更多的关注非理想采集虹膜的识别。非理想采集虹膜图像由于存在反光、睫毛和眼睑遮挡等干扰，图像质量大为降低，这为虹膜边界的准确定位带来了很大困难。为此，人们提出了不同的方法。例如，等利用k均值模糊聚类实现对虹膜图像灰度的特征分类；Pundlik等利用图割实现虹膜分割，但是图割方法依赖于像素间的统计关系，忽略了边界信息，且迭代过程较长；He等提出了一种混合分割方法，即首先利用线性插值算法消除反光，再由弦长均衡方法搜索瞳孔中心；Jarjes等利用snake模型和角积分投影方法实现虹膜分割；Roy通过水平集方法实现虹膜边界的提取。然而，由于Snake和水平集两种方法均依赖于由偏微分方程控制的曲线演化模型，此类虹膜分割方法收敛速度较慢，且对曲线的初始化位置非常敏感。

发明内容

本发明针对上述关于非理想采集虹膜分割方法存在的不足，提出了一种基于非线性尺度空间的虹膜分割算法。这一算法分为两个阶段，能够在消除反光和睫毛遮挡的基础上，分别实现内边界定位和外边界定位。

本发明采用的算法方案如下：

一种基于非线性尺度空间的虹膜图像分割算法，包括以下步骤：

步骤1对虹膜内边界定位：利用非线性尺度空间所具有的尺度演化特性，对原始虹膜图像进行平滑处理，消除瞳孔区域的亮斑，保证瞳孔区域的灰度处于整幅图像的最低水平，再通过阈值实现瞳孔检测，实现虹膜内边界定位；

步骤2对虹膜外边界定位：利用非线性尺度空间的尺度演化特性，对原始虹膜图像进行平滑处理，消除几何尺度较小的睫毛对虹膜区域的遮挡，实现了对外边界的准确定位。

所述的步骤1具体包括以下步骤：

（1-1）对输入的虹膜图像I进行框架Laplace金字塔分解，得到概貌图像Ic和细节图像Id；

（1-2）对概貌图像Ic进行非线性尺度空间下的各向异性扩散，得到图像Ic_diffu1；

（1-3）考虑到正常情况下瞳孔区域的灰度值最低，则设定阈值Th，将图像Ic_diffu1二值化，低于此阈值的像素的灰度值置为0，得到L区域；否则置1，得到H区域；

（1-4）对L区域进行闭运算，得到图像Ic_pupil，即概貌图像Ic的二值化的瞳孔区域，其中闭运算采用的结构算子为圆，用b表示，即

$\mathrm{Ic}\_\mathrm{pupil}=(L\⊕b)\mathrm{\Θb}---\left(1\right)$

（1-5）利用插值算法将图像Ic_pupil映射为与虹膜图像I同样大小的图像，得到图像I_pupil，采用双线性插值实现；

（1-6）以图像I_pupil的L区域为模板，覆盖原始虹膜图像I中对应的区域，得到消除反光后的图像I_r；

（1-7）计算图像I_r瞳孔区域的质心C，即

$C=(C_x,C_y)=(\frac{1}{N_{\mathrm{pupil}}}\underset{I\_r(i,j)\∈L}{\mathrm{\Σ}}i,\frac{1}{N_{\mathrm{pupil}}}\underset{I\_r(i,j)\∈L}{\mathrm{\Σ}}j)---\left(2\right)$

其中(i,j)为图像像素点坐标，N_pupil为瞳孔区域的像素数量；又以C点为圆心，对I_r进行极坐标变换和边缘检测，以半径方向上幅值响应最大的点为内边界点；统计内边界点到质心C的距离，将出现概率最大的距离Ri作为内边界的半径。

所述的步骤（1-1）所述的对输入的虹膜图像I进行框架Laplace金字塔分解，这一过程以矩阵形式表示，其一层分解公式为

$\left(\begin{array}{c}\mathrm{Ic}\\ \mathrm{Id}\end{array}\right)=\left(\begin{array}{c}(\↓2)H\\ I-G(\↑2)(\↓2)H\end{array}\right)I---\left(3\right)$

其中H和G分别表示满足完全重构条件(即框架结构)的低通和高通滤波器，Ic表示概貌图像，Id表示细节图像；

所述的步骤（1-2）所述的异性扩散，采用的扩散控制函数选择（4）式

$g\left(\right|\▿I\left|\right)=\frac{1}{1+{\left(\frac{|\▿I|}{K}\right)}^2}---\left(4\right)$

是扩散控制函数，它关于梯度幅度单调递减，K为设定的常数，梯度表示图像边缘的检测函数。

得到图像Ic_diffu1，其数值实现由(5)式确定，参数设定为λⁿ＝0.25；

$I^{n+1}(i,j)=I^n(i,j)+{\mathrm{\λ}}^n\·g\left(\right|\▿I^n(i,j)\left|\right)(I_{\mathrm{xx}}^n(i,j)+I_{\mathrm{yy}}^n(i,j))$

$+{\mathrm{\λ}}^n\·g^{\′}\left(\right|\▿I^n(i,j)\left|\right)\left(\frac{[I_x^n(i,j){]}^2\·I_{\mathrm{xx}}^n(i,j)+[I_y^n(i,j){]}^2\·I_{\mathrm{yy}}^n(i,j)}{|\▿I^n(i,j){|}^2}\right)---\left(5\right)$

所述的步骤（1-3）用如下公式表示：

$\mathrm{Ic}\_b(i,j)=\left\{\begin{array}{c}H,\mathrm{if\; Ic}\_\mathrm{difful}(i,j)\&GreaterEqual;\mathrm{Th}\\ L,\mathrm{if\; Ic}\_\mathrm{difful}(i,j)<\mathrm{Th}\end{array}\right.---\left(6\right)$

所述的步骤（1-4）主要消除小尺度的反光区域；二值化操作的目的是得到瞳孔区域，其中瞳孔区域为黑色，背景及反光区域为白色，再将此区域进行形态学闭运算，即膨胀+腐蚀，瞳孔内部的反光被黑色填充，即去反光操作，得到的L区域(Ic_pupil)即为概貌图像中的瞳孔区域。

所述的步骤（1-7）的最大的距离Ri表示为：

$\mathrm{Ri}=\underset{k\&Element;K}{\arg \max }{p}_k---\left(7\right)$

其中p_k＝p(r＝Ri_k)，表示半径为Ri_k的概率，K为所有存在的半径数量；Ri_k∈RI,RI为所有检测出的半径长度的集合。

所述的步骤2对虹膜外边界定位，包括以下步骤；

（2-1）将由步骤1得到的消除反光后的图像I_r进行Laplace金字塔分解，得到概貌图像I_rc，内边界圆心C在概貌图像I_rc中被映射为点分解过程与(7)式相同；

I(x,y,t)＝I₀(x,y)*G(x,y,t)(8)

其中I₀(x,y)和I(x,y,t)表示初始时刻和t时刻的图像，G(x,y,t)为t时刻的高斯核函数，其方差为时间变量t；

（2-2）对概貌图像I_rc进行非线性扩散，得到图像I_rc_diffu2，参数设定为Δt＝0.25和λⁿ＝0.25；

（2-3）判断与图像I_rc_diffu2左右边界的距离，若与右边界距离较大，则以为圆心、在角度[π,7π/6]范围内进行极坐标变换；否则以为圆心、在[11π/6,2π]范围内进行极坐标变换；此处变换所取半径为与左右边界距离的较大者，由此步骤得到的(极坐标变换)图像的瞳孔以外的部分记为图像I_rp；

（2-4）将图像I_rp按半径方向等分三部分，对距瞳孔最近的部分(用I_rp_upper表示)进行极大值滤波，最远的部分(用I_rp_bottom表示)进行极小值滤波，中间部分(用I_rp_middle表示)进行中值滤波，得到图像

（2-5）在图像中，基于区域平均灰度差的曲线演化策略，以此实现外边界的定位。

(2-6)将上述迭代得到的边界C₀映射到图像I_r，即为虹膜外边界。

所述的步骤（2-2）所述的扩散控制函数选择(4)式，

$g\left(\right|\&dtri;I\left|\right)=\frac{1}{1+{\left(\frac{|\&dtri;I|}{K}\right)}^2}---\left(4\right)$

所述的步骤（2-5）的迭代过程表示为

a)将中值滤波区域中靠近最小值滤波区域的边界设定为初始外边界，，记为C₀，同时得到由C₀划分的两个区域A_i和A_e；所述的A_i为靠近瞳孔区域，A_e为远离瞳孔区域；

b)计算由C₀划分的两个区域的平均灰度，其中靠近瞳孔的记为d₀，远离瞳孔的记为d₀'，即：

$\left\{\begin{array}{c}{d}_0=\mathrm{mean}\left(A_i\right)\\ d_0^{\&prime;}=\mathrm{mean}\left(A_e\right)\end{array}\right.$

c)将C₀向远离瞳孔的方向移动k个像素，得到新的边界C₁以及区域A_i ¹和A_e ¹，计算以C₁划分的区域的新的平均灰度d和d'；这一过程表示为：

$\left\{\begin{array}{c}d=\mathrm{mean}\left({A_i}^1\right)\\ d^{\&prime;}=\mathrm{mean}\left({A_e}^1\right)\end{array}\right.$

d)按照|d₀-d'₀|与|d-d'|的大小关系确定边界是否继续移动，即若|d₀-d'₀|≤|d-d'|，令C₀＝C₁,d₀＝d,d'₀＝d'并返回步骤b)，否则输出C₀作为外边界定位结果。

本发明的有益效果如下：

(1)本算法将非线性尺度空间引入虹膜预处理，实现了对虹膜内边界(即瞳孔边界)的准确定位。虹膜内边界定位的难点在于瞳孔的检测。瞳孔在虹膜图像中的灰度水平最低，一般采取阈值策略进行检测。本算法利用非线性尺度空间所具有的尺度演化特性，对原始虹膜图像进行平滑，消除了瞳孔区域的亮斑，保证了瞳孔区域的灰度处于整幅图像的最低水平，再通过阈值实现瞳孔检测。

(2)本算法在外边界定位过程中，同样利用非线性尺度空间的尺度演化特，对原始虹膜图像进行平滑，消除了几何尺度较小的睫毛对虹膜区域的遮挡，实现了对外边界的准确定位。

附图说明

图1虹膜内边界定位流程图；

图2虹膜外边界定位流程图；

图3非线性尺度空间对原始虹膜图像进行平滑的结果；

图4本算法处理结果一；

图5本算法处理结果二。

具体实施方式

下面结合附图对本发明进行详细说明：

首先针对本发明所采用的非线性尺度空间算法进行介绍，非线性尺度空间的概念源自线性尺度空间,也称为高斯尺度空间，其演化方程可以表述为各向同性扩散(IsotropicDiffusion)，即

I(x,y,t)＝I₀(x,y)*G(x,y,t)(9)

其中I₀(x,y)和I(x,y,t)表示初始时刻和t时刻的图像，G(x,y,t)为t时刻的高斯核函数，其方差为时间变量t。此外，线性尺度空间还有另一种等价形式，即

$\left\{\begin{array}{c}\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(\&dtri;I)=\mathrm{\&Delta;I}(x,y)\\ I(x,y,0)=I_0(x,y)\end{array}---\left(10\right)\right.$

这里Δ为Laplace算子，div(·)表示散度。由(1)式可以看出，线性尺度空间描述的图像关于时间的演化规律是：初始图像被参数为t的二维高斯核函数平滑，其包含的细节信息随时间的增加而不断减少。因此，这一方法的最大不足是造成图像细节的损失。

与线性尺度空间相比，非线性尺度空间下的各向异性扩散(AnisotropicDiffusion)在最大优势在于平滑图像的同时能够有效保护边缘。因此这一方法在图像处理和计算机视觉领域得到了广泛应用。Perona和Malik^[14]最早提出了基本形式的非线性扩散方程，即P-M模型。作为对(2)式的泛化形式，P-M模型表示为如下的偏微分方程

$\frac{\&PartialD;I}{\&PartialD;t}=\mathrm{div}(g\left(\right|\&dtri;I\left|\right)\&dtri;I)=g\left(\right|\&dtri;I\left|\right)\&CenterDot;\mathrm{\&Delta;I}+\&lang;\&dtri;g,\&dtri;I\&rang;---\left(11\right)$

其中I为t时刻的图像，<·>表示内积运算。是扩散控制函数，它关于梯度幅度单调递减，其表达式为

$g\left(\right|\&dtri;I\left|\right)=\exp (-{\left(\frac{|\&dtri;I|}{K}\right)}^2)---\left(12\right)$

或

$g\left(\right|\&dtri;I|=)\frac{1}{1+{\left(\frac{|\&dtri;I|}{K}\right)}^2}---\left(4\right)$

K为设定的常数。就非线性扩散而言，梯度实际上是图像边缘的检测函数，而则被设为关于梯度的单调递减(或非增)函数。P-M方程对图像的平滑体现了非线性扩散的基本思想：越大，表明当前像素点处存在强边缘，故其平滑程度越小；反之，则表明当前像素点处于平坦区域，平滑程度越大。换言之，非线性扩散具有良好的边缘保护特性。另外，非线性扩散的一个重要性质是其尺度演化(平滑)特性：随着演化即时间的推进，具有小尺度特征的几何或纹理结构会被逐渐消除，而大尺度的边缘信息得以保留^[15]。随着细小纹理结构的消失，图像内容被大大简化，加之主要边界得以保护，这为图像分割提供了很大帮助。图像的非线性演化过程由迭代实现。若函数(a为常数)，则P-M方程退化为线性的各向同性扩散^[14]。

对于二维数字图像，由(3)式及梯度下降法原理，可得P-M模型的(4邻域)数值实现为 $I^{n+1}(i,j)=I^n(i,j)+{\mathrm{\&lambda;}}^n\&CenterDot;g\left(\right|\&dtri;I^n(i,j)\left|\right)(I_{\mathrm{xx}}^n(i,j)+I_{\mathrm{yy}}^n(i,j))$

$+{\mathrm{\&lambda;}}^n\&CenterDot;g^{\&prime;}\left(\right|\&dtri;I^n(i,j)\left|\right)\left(\frac{[I_x^n(i,j){]}^2\&CenterDot;I_{\mathrm{xx}}^n(i,j)+[I_y^n(i,j){]}^2\&CenterDot;I_{\mathrm{yy}}^n(i,j)}{|\&dtri;I^n(i,j){|}^2}\right)---\left(5\right)$

其中 $|\&dtri;I^n(i,j)|=\sqrt{[I_x^n(i,j){]}^2+[I_y^n(i,j){]}^2}$ 并且

$\left\{\begin{array}{c}{I}_x^n(i,j)=I^n(i+1,j)-I^n(i,j)\\ I_y^n(i,j)=I^n(i,j+1)-I^n(i,j)\\ I_{\mathrm{xx}}^n(i,j)=I^n(i+1,j)+I^n(i-1,j)-2I^n(i,j)\\ I_{\mathrm{yy}}^n(i,j)=I^n(i,j+1)+I^n(i,j-1)-2I^n(i,j)\end{array}\right.$

这里Iⁿ(i,j)表示第n次迭代后图像在(i,j)处的像素值，和分别表示图像I在(i,j)处水平和垂直方向的一阶差分，和分别表示二阶差分。此外λⁿ为第n次迭代时的权值，迭代次数可以根据实际需要人为设定。图-1展示了虹膜图像在P-M模型下的演化结果，参数设定为λⁿ≡0.25，扩散函数g(·)取(5)式。很明显，小尺度的睫毛和眼睑随着迭代次数的增加而愈加模糊，而大尺度的瞳孔则得以保留。

本发明中基于非线性尺度空间的虹膜图像分割算法，主要过程如下：

第一阶段：内边界定位

虹膜内边界的提取是整个边界定位算法的关键。就虹膜图像而言，算法必须能够克服各种常见的尤其是睫毛的干扰，保证内边界和瞳孔中心点的准确定位，同时在保证算法可靠性的基础上尽量节省计算。由此本文提出内虹膜图像边界算法由图-1所示，具体步骤由如下：

(1-1)对输入的虹膜图像I进行框架Laplace金字塔分解^[16]，这一过程以矩阵形式表示，其一层分解公式为

$\left(\begin{array}{c}\mathrm{Ic}\\ \mathrm{Id}\end{array}\right)=\left(\begin{array}{c}(\&DownArrow;2)H\\ I-G(\&UpArrow;2)(\&DownArrow;2)H\end{array}\right)I---\left(3\right)$

其中H和G分别表示满足完全重构条件(即框架结构)的低通和高通滤波器，得到概貌图像Ic和细节图像Id；

(1-2)对Ic进行非线性尺度空间下的各向异性扩散，扩散控制函数选择(11)式，得到图像Ic_diffu1，其数值实现由(15)式确定，参数设定为λⁿ=0.25；

(1-3)考虑到一般情况下瞳孔区域的灰度值最低，则设定阈值Th，将Ic_diffu1(Ic_diffu1为步骤(1-2)的结果)二值化，低于此阈值的像素的灰度值置为0，得到L区域；否则置1，得到H区域，即得到二值图像Ic_b

$\mathrm{Ic}\_b(i,j)=\left\{\begin{array}{c}H,\mathrm{if\; Ic}\_\mathrm{difful}(i,j)\&GreaterEqual;\mathrm{Th}\\ L,\mathrm{if\; Ic}\_\mathrm{difful}(i,j)<\mathrm{Th}\end{array}\right.---\left(6\right)$

(i,j)表示图像坐标。

(1-4)对L区域进行闭运算，得到图像Ic_pupil，即概貌图像的二值化的瞳孔区域，其中闭运算采用的结构算子为圆，用b表示，即

$\mathrm{Ic}\_\mathrm{pupil}=(L\&CirclePlus;b)\mathrm{\&Theta;b}---\left(1\right)$

这一运算将消除小尺度的反光区域。二值化操作的目的是得到瞳孔区域(瞳孔区域为黑色，背景及反光区域为白色)，再将此区域进行形态学闭运算(膨胀+腐蚀)，瞳孔内部的反光被黑色填充，即去反光操作，得到的L区域(Ic_pupil)即为概貌图像中的瞳孔区域；

(1-5)利用插值算法将Ic_pupil映射为与I同样大小的图像，得到I_pupil，采用双线性插值实现；

(1-6)以I_pupil的L区域为模板，覆盖原始图像I中对应的区域，得到消除反光后的图像I_r；

(1-7)计算I_r瞳孔区域的质心C，即

$C=(C_x,C_y)=(\frac{1}{N_{\mathrm{pupil}}}\underset{I\_r(i,j)\&Element;L}{\mathrm{\&Sigma;}}i,\frac{1}{N_{\mathrm{pupil}}}\underset{I\_r(i,j)\&Element;L}{\mathrm{\&Sigma;}}j)---\left(2\right)$

其中N_pupil为瞳孔区域的像素数量。又以C点为圆心，对I_r进行极坐标变换和边缘检测，以半径方向上幅值响应最大的点为内边界点；统计内边界点到质心C的距离，将出现概率最大的距离Ri作为内边界的半径。这一过程可表示为

$\mathrm{Ri}=\underset{k\&Element;K}{\arg \max }{p}_k---\left(7\right)$

其中p_k＝p(r＝Ri_k)，表示半径为Ri_k的概率；Ri_k∈RI,RI为所有检测出的半径长度的集合。由上述过程得到的质心C即为瞳孔中心，亦即虹膜内边界的圆心；Ri为内边界的半径。

第二阶段：外边界定位

虹膜外边界通常表现为变化平缓的弱边缘，故常用的边缘检测算子对其检测的效果并不理想。本文采用基于区域统计信息(平均灰度差)的轮廓演化方法实现外边界的定位(见图2)，具体算法如下：

(2-1)将由步骤(1-6)得到的消除反光后的图像I_r进行Laplace金字塔分解，得到概貌图像I_rc，内边界圆心C在I_rc中被映射为点分解过程与(7)式相同；

(2-2)对I_rc进行非线性扩散，得到图像I_rc_diffu2，扩散控制函数选择(11)式(见附录)，参数设定为Δt＝0.25和λⁿ＝0.25；

(2-3)判断与图像I_rc_diffu2左右边界的距离，若与右边界距离较大，则以为圆心、在角度[π,7π/6]范围内进行极坐标变换；否则以为圆心、在[11π/6,2π]范围内进行极坐标变换。此处变换所取半径为与左右边界距离的较大者。由此步骤得到的(极坐标变换)图像的瞳孔以外的部分记为图像I_rp；

(2-4)将图像I_rp按半径方向等分三部分，对距瞳孔最近的部分(用I_rp_upper表示)进行极大值滤波，最远的部分(用I_rp_bottom表示)进行极小值滤波，中间部分(用I_rp_middle表示)进行中值滤波，得到图像这一过程可以表示为

$\stackrel{\&OverBar;}{I\_\mathrm{rp}}(i,j)=\left\{\begin{array}{c}\min (I\_\mathrm{rp}(i,j)),(i,j)\&Element;I\_r{p}_{\mathrm{upper}}\\ \mathrm{med}(I\_\mathrm{rp}(i,j)),(i,j)\&Element;I\_r{p}_{\mathrm{middle}}\\ \max (I\_\mathrm{rp}(i,j)),(i,j)\&Element;I\_r{p}_{\mathrm{bottom}}\end{array}\right.---\left(13\right)$

并且(i,j)满足(i,j)∈Ω_(i,j)，即(i,j)属于以自身为中心的邻域。外边界定位由基于区域信息的曲线演化(迭代)实现，这一过程可以看作简化的几何活动轮廓模型。为防止边界演化陷入局部极值，我们需要增强虹膜与巩膜之间的对比度同时不破坏二者之间的边界，为此采用上述分段序统计滤波的策略。这一过程见图2；

(2-5)在图像中，我们提出了基于区域平均灰度差的曲线演化策略，以此实现外边界的定位。迭代过程可表示为：

a)将中值滤波区域中靠近最小值滤波区域的边界设定为初始外边界，，记为C₀，同时得到由C₀划分的两个区域A_i(靠近瞳孔区域)和A_e(远离瞳孔区域)；

b)计算由C₀划分的两个区域的平均灰度，其中靠近瞳孔的记为d₀，远离瞳孔的记为d₀'，即

$\left\{\begin{array}{c}{d}_0=\mathrm{mean}\left(A_i\right)\\ d_0^{\&prime;}=\mathrm{mean}\left(A_e\right)\end{array}\right.$

c)将C₀向远离瞳孔的方向移动k个像素，得到新的边界C₁以及区域A_i ¹和A_e ¹，计算以C₁划分的区域的新的平均灰度d和d'。这一过程表示为

$\left\{\begin{array}{c}d=\mathrm{mean}\left({A_i}^1\right)\\ d^{\&prime;}=\mathrm{mean}\left({A_e}^1\right)\end{array}\right.$

d)按照|d₀-d'₀|与|d-d'|的大小关系确定边界是否继续移动，即若|d₀-d'₀|≤|d-d'|，令C₀＝C₁,d₀＝d,d'₀＝d'并返回步骤b)，否则输出C₀作为外边界定位结果；

(S2-6)将上述迭代得到的边界C₀映射到图像I_r，即为虹膜外边界。

如图3所示为非线性尺度空间对原始虹膜图像进行平滑的结果，从左向右平滑尺度逐渐增加；如图4所示为本算法处理结果一，左列为原始图像，中列为内边界定位结果，右列为外边界定位结果；如图5所示为本算法处理结果二，左列为原始图像，中列为内边界定位结果，右列为外边界定位结果。在近红外波段采集的仅包含眼部区域的虹膜图像，以通用格式(如bmp,jpg,tiff等)存储在计算机内，调用本算法进行处理，再将处理结果用识图软件显示即可。

Claims (8)

1.一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，包括以下步骤：

步骤1对虹膜内边界定位：利用非线性尺度空间所具有的尺度演化特性，对原始虹膜图像进行平滑处理，消除瞳孔区域的亮斑，保证瞳孔区域的灰度处于整幅图像的最低水平，再通过阈值实现瞳孔检测，实现虹膜内边界定位；

步骤2对虹膜外边界定位：利用非线性尺度空间的尺度演化特性，对原始虹膜图像进行平滑处理，消除几何尺度较小的睫毛对虹膜区域的遮挡，实现了对外边界的准确定位；

所述的步骤1具体包括以下步骤：

(1-1)对输入的虹膜图像I进行框架Laplace金字塔分解，得到概貌图像Ic和细节图像Id；

(1-2)对概貌图像Ic进行非线性尺度空间下的各向异性扩散，得到图像Ic_diffu1；

(1-3)考虑到正常情况下瞳孔区域的灰度值最低，则设定阈值Th，将图像Ic_diffu1二值化，低于此阈值的像素的灰度值置为0，得到L区域；否则置1，得到H区域；

(1-4)对L区域进行闭运算，得到图像Ic_pupil，即概貌图像Ic的二值化的瞳孔区域，其中闭运算采用的结构算子为圆，用b表示，即

$Ic\_pupil=(L\&CirclePlus;b)\mathrm{\&Theta;}b---\left(1\right)$

(1-5)利用插值算法将图像Ic_pupil映射为与虹膜图像I同样大小的图像，得到图像I_pupil，采用双线性插值实现；

(1-6)以图像I_pupil的L区域为模板，覆盖原始虹膜图像I中对应的区域，得到消除反光后的图像I_r；

(1-7)计算图像I_r瞳孔区域的质心C，即

$C=(C_x,C_y)=(\frac{1}{N_{pupil}}\underset{I\_r(i,j)\&Element;L}{\&Sigma;}i,\frac{1}{N_{pupil}}\underset{I\_r(i,j)\&Element;L}{\&Sigma;}j)---\left(2\right)$

其中(i,j)为图像像素点坐标，N_pupil为瞳孔区域的像素数量；又以C点为圆心，对I_r进行极坐标变换和边缘检测，以半径方向上幅值响应最大的点为内边界点；统计内边界点到质心C的距离，将出现概率最大的距离Ri作为内边界的半径。

2.如权利要求1所述的一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，所述的步骤(1-2)所述的异性扩散，采用的扩散控制函数选择(4)式

$g\left(\right|\&dtri;I\left|\right)=\frac{1}{1+{\left(\frac{|\&dtri;I|}{K}\right)}^2}---\left(4\right)$

g(|▽I|)是扩散控制函数，它关于梯度幅度|▽I|单调递减，K为设定的常数，梯度|▽I|表示图像边缘的检测函数；

得到图像Ic_diffu1，其数值实现由(5)式确定，参数设定为λⁿ＝0.25；

$\begin{array}{c}{I}^{n+1}(i,j)=I^n(i,j)+{\mathrm{\&lambda;}}^n\&CenterDot;g\left(|\&dtri;I^n\left(i,j\right)|\right)(I_{xx}^n(i,j)+I_{yy}^n(i,j))\\ +{\mathrm{\&lambda;}}^n\&CenterDot;g^{\&prime;}\left(|\&dtri;I^n\left(i,j\right)|\right)\left(\frac{{\&lsqb;I_x^n(i,j)\&rsqb;}^2\&CenterDot;I_{xx}^n(i,j)+{\&lsqb;I_y^n(i,j)\&rsqb;}^2\&CenterDot;I_y^n(i,j)}{|\&dtri;I^n\left(i,j\right){|}^2}\right)\end{array}---\left(5\right)$

其中I(i,j)表示图像在坐标(i,j)处的灰度值，(i,j)表示二维坐标，Iⁿ(i,j)表示第n次迭代后图像的灰度值，Iⁿ⁺¹(i,j)表示n+1次迭代后图像的灰度值，▽Iⁿ(i,j)表示第n次迭代后图像的梯度，|▽Iⁿ(i,j)|表示▽Iⁿ(i,j)的绝对值，g(|▽Iⁿ(i,j)|)表示以▽Iⁿ(i,j)为自变量的单调非增函数，g'(|▽Iⁿ(i,j)|)表示g(|▽Iⁿ(i,j)|)的导数，表示第n次迭代后图像沿水平方向的二阶偏导数，表示第n次迭代后图像沿垂直方向的二阶偏导数，表示第n次迭代后图像沿水平方向的一阶偏导数，表示第n次迭代后图像沿垂直方向的一阶偏导数。

3.如权利要求1所述的一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，所述的步骤(1-3)用如下公式表示：

$Ic\_b(i,j)=\left\{\begin{array}{ccc}H,& if& Ic\_difu1(i,j)\&GreaterEqual;Th\\ L,& if& Ic\_difu1(i,j)<Th\end{array}\right.---\left(6\right)$

Ic_b(i,j)表示对图像Ic_diffu1进行如(6)式运算后的到的图像。

4.如权利要求1所述的一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，所述的步骤(1-4)主要消除小尺度的反光区域；二值化操作的目的是得到瞳孔区域，其中瞳孔区域为黑色，背景及反光区域为白色，再将此区域进行形态学闭运算，即膨胀+腐蚀，瞳孔内部的反光被黑色填充，即去反光操作，得到的L区域即为概貌图像中的瞳孔区域。

5.如权利要求1所述的一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，所述的步骤(1-7)的最大的距离Ri表示为：

$Ri=\underset{k\&Element;K}{\mathrm{argmax}}{p}_k---\left(7\right)$

其中p_k＝p(r＝Ri_k)，表示半径为Ri_k的概率，K为所有存在的半径数量；Ri_k∈RI,RI为所有检测出的半径长度的集合。

6.如权利要求5所述的一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，所述的步骤2对虹膜外边界定位，包括以下步骤；

(2-1)将由步骤1得到的消除反光后的图像I_r进行Laplace金字塔分解，得到概貌图像I_rc，内边界圆心C在概貌图像I_rc中被映射为点

(2-2)对概貌图像I_rc进行非线性扩散，得到图像I_rc_diffu2；

(2-3)判断与图像I_rc_diffu2左右边界的距离，若与右边界距离较大，则以为圆心、在角度[π,7π/6]范围内进行极坐标变换；否则以为圆心、在[11π/6,2π]范围内进行极坐标变换；此处变换所取半径为与左右边界距离的较大者，由此步骤得到的图像的瞳孔以外的部分记为图像I_rp；

(2-4)将图像I_rp按半径方向等分三部分，对距瞳孔最近的部分进行极大值滤波，最远的部分进行极小值滤波，中间部分进行中值滤波，得到图像

(2-5)在图像中，基于区域平均灰度差的曲线演化策略，以此实现外边界的定位；

(2-6)将步骤(2-5)得到的边界C₀映射到图像I_r，即为虹膜外边界。

7.如权利要求6所述的一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，所述的步骤(2-2)对概貌图像I_rc进行非线性扩散，其扩散控制函数采用公式(4)，

$g\left(\right|\&dtri;I\left|\right)=\frac{1}{1+{\left(\frac{|\&dtri;I|}{K}\right)}^2}---\left(4\right)$

其中|▽I|表示图像梯度的绝对值，K表示扩散因子，g(|▽I|)表示以|▽I|为自变量的单调非增函数。

8.如权利要求6所述的一种基于非线性尺度空间的虹膜图像分割方法，其特征在于，所述的步骤(2-5)表示为：

a)将中值滤波区域中靠近最小值滤波区域的边界设定为初始外边界，记为C₀，同时得到由C₀划分的两个区域A_i和A_e；所述的A_i为靠近瞳孔区域，A_e为远离瞳孔区域；

b)计算由C₀划分的两个区域的平均灰度，其中靠近瞳孔的记为d₀，远离瞳孔的记为d₀'，即：

$\left\{\begin{array}{c}{d}_0=mean\left(A_i\right)\\ d_0^{\&prime;}=mean\left(A_e\right)\end{array}\right.$

c)将C₀向远离瞳孔的方向移动k个像素，得到新的边界C₁以及区域A_i ¹和A_e ¹，计算以C₁划分的区域的新的平均灰度d和d'；这一过程表示为：

$\left\{\begin{array}{c}d=mean\left({A_i}^1\right)\\ d^{\&prime;}=mean\left({A_e}^1\right)\end{array}\right.$

d)按照|d₀-d'₀|与|d-d'|的大小关系确定边界是否继续移动，即若|d₀-d'₀|≤|d-d'|，令C₀＝C₁,d₀＝d,d'₀＝d'并返回步骤b)，否则输出C₀作为外边界定位结果。