一种基于混合高斯模型的虹膜图像分割算法

摘要

本发明公开了一种基于混合高斯模型的虹膜图像分割算法，包括内边界定位、外边界定位和眼睑检测的步骤；本发明本算法将混合高斯模型引入虹膜预处理，实现了对虹膜内边界的准确定位。虹膜内边界定位的难点在于瞳孔的检测。瞳孔在虹膜图像中的灰度水平最低，一般采取阈值策略进行检测。本算法利用混合高斯模型所具有的对不同数据的自适应统计特性，对原始虹膜图像进行分割，获得灰度最低的瞳孔区域；然后在瞳孔内部初始化虹膜中心，利用弦长搜索策略实现内边界定位。本算法在外边界定位过程中，同样在混合高斯模型分割的基础上，通过使局部对比度最大化实现了对外边界的准确定位。

说明书

技术领域

本发明涉及一种基于混合高斯模型的虹膜图像分割算法。 

背景技术

生物特征识别(Biometrics)是指通过信息技术，利用人类自身的生理或行为特征实现个人身份认定。这些特征包括虹膜、指纹、掌纹、人脸、声音等人体的生理特征，以及签名动作、行走步态、敲击键盘的力度等行为特征。由于这些生理或行为特征因人而异且具有相当的稳定性，因此成为鉴别身份的可靠依据。近年来，生物特征识别技术已在信息安全、金融交易、社会安全、人员管理、医疗卫生等领域获得了广泛应用[1]。 

作为一项重要的生物特征识别技术，虹膜识别得到了科学界和工业界的高度关注[2,3]，原因在于与其他生物特相比，虹膜唯一性高、精度高、稳定性好，适用人群广泛。目前，一个典型的虹膜识别系统主要由虹膜图像预处理、特征提取、归一化和编码几部分组成。其中虹膜图像预处理的主要任务是正确分割图像中的虹膜区域，即定位虹膜内外边界、消除睫毛和眼睑遮挡，从而为后续的特征提取和编码过程提供有效信息。预处理是整个识别过程的基础和关键，它会对最终的识别效果产生重要影响。最早的虹膜边界定位方法由Daugman提出，即经典的积分/微分算子(Integro-differential operator)^[4]；Wildes提出的算法将边缘检测与Hough变换相结合^[5]；Boles等采用一维三次样条小波提取出的图像过零点作为特征^[6]。随后，人们在虹膜图像预处理阶段引入了多种工具，如活动轮廓模型、傅里叶变换与几何投影、统计方法、Gabor滤波器等^[7,8]，取得了良好效果。近年来，考虑到虹膜识别方法的实用性，人们开始更多的关注虹膜的识别^{[9,10,11,12,13]}。虹膜图像由于存在反光、睫毛和眼睑遮挡等干扰，图像质量大为降低，这为虹膜边界的准确定位带来了很大困难。因此，虹膜识别的关键在于虹膜的内外边界定位。 

考虑到虹膜识别系统的实用性，研究者已经对虹膜预处理给予很大关注。由于存在反光、睫毛和眼睑遮挡等干扰，虹膜边界的准确定位存在很大困难。为此，人们近年来提出了许多方法。例如，Proenca^[9]等提出利用k均值模糊聚类对虹膜图像灰度的特征分类；Pundlik等^[10]通过图割分割虹膜，缺点是图割依赖于像素间的统计关系，忽略了边界信息，且迭代过程较长；He等^[11]提出了一种混合分割方法，即首先利用线性插值技术消除反光，再由弦长均衡方法搜索瞳孔中心；Jarjes等^[12]利用snake模型和角积分投影进行虹膜分割；Roy等^[13]通过 水平集演化实现虹膜边界提取。需要指出，由于Snake和水平集方法均依赖于由偏微分方程控制的曲线演化模型，此类方法收敛速度慢，且对曲线的初始化位置非常敏感。此外，Verma等^[14]提出利用图像增强技术分割虹膜。 

参考文献： 

[1]A K Jain,A Ross,S Pankanti.Biometrics:a tool for information security[J].IEEE Trans.Information Forensics and Security,2006,1(2),125-143. 

[2]J Daugman.How iris recognition works[J].IEEE Trans.Circuits and Systems for Video Technology,2004,14(1):21-30. 

[3]K W Bowyer,K Hollingsworth,P J Flynn.Image-understanding for iris biometrics:a survey[J].Computer Vision and Image Understanding,2008,110(2):281-307. 

[4]J Daugman.High confidence visual recognition of persons by a test of statistical independence[J].IEEE Trans.Pattern Analysis and Machine Intelligence,1993,15(11):1148-1161. 

[5]R P Wildes.Iris recognition:an emerging biometrics technology[A].Proc.IEEE[C],1997,85(9):1348-1363. 

[6]W Boles,B Boashah.A human identification technique using images of the iris and wavelet transform[J].IEEE Trans.Signal Processing,1998,46(4):1185-1188. 

[7]J Daugman.New methods in iris recognition[J].IEEE Trans.System,Man,and Sybernetics—Part B:Cybernetics,2007,37(5):1167-1175. 

[8]L Yu,D Zhang,K Wang.The relative distance of key point based iris recognition[J].Pattern Recognition,2007,40(2):423-430. 

[9]HL A Alexandre.Iris segmentation methodology for non-cooperative recognition [A].IEE Proc.VISP[C],2006,153(2):199-205. 

[10]S J Pundlik,D L Woodard,S T Birchfield.Non-ideal iris segmentation using graph cuts [A].Proc.IEEE CVPR[C],2008:1-6. 

[11]Z He,T Tan,Z Sun,X Qiu.Toward accurate and fast iris segmentation for iris biometrics [J].IEEE Trans.Pattern Analysis and Machine Intelligence,2009,31(9):1670-1684. 

[12]A A Jarjes,K Wang,G J Mohammed.A new iris segmentation method based on improved snake model and angular integral projection.Research Journal of Applied Science,Engineering and Technology,2011,3(6):558-568. 

[13]K Roy,P Bhattacharya,C Y Suen.Iris segmentation using variational level set method.Optics and Laser in Engineering,2011,49(1):578-588. 

[14]A Verma,C Liu,J Jia.Iris recognition based on robust iris segmentation and image enhancement,International Journal of Biometrics,2012,4(1):56-75. 

[15]N Shental,A Bar-Hillel,T Hertz,D Weinshall.Computing Gaussian mixture models with EM using equivalence constraints.Advances in Neural Information Processing Systems16,2004,16(8),465-472. 

[16]A Dempster,N Laird,D Robin.Maximum likelihood from incomplete data via the EM algorithm.Journal of the Royal Statistical Society,Series B,1977,39(1):1-38. 

发明内容

本算法提出一种基于混合高斯模型^[15](Gaussian Mixture Model,GMM)的虹膜分割方法。混合高斯模型通过多个高斯概率概率函数之和，在某种准则下实现对图像直方图的最佳逼近。由于这一方法以图像自身的统计特征(灰度的概率分布)为逼近对象，对不同图像具有很强的适应性，因此广泛应用于图像分割。在本算法方法中，在内边界定位阶段我们首先利用混合高斯模型实现对虹膜图像的建模并对其进行分割，从而实现对瞳孔的检测。考虑到瞳孔形状经常因为睫毛和眼睑遮挡而不规则，我们由弦长均衡策搜索其中心(亦即虹膜中心)坐标并确定内圆(内边界)半径。在外边界定位过程中，我们根据虹膜外边界灰度变化缓慢的特点，提出了基于区域灰度差的曲线演化策略，通过迭代搜索外边界。为确保曲线能够全局收敛，我们利用序统计滤波虹膜进行处理。最后，我们利用二次曲线实现眼睑检测。 

为了实现上述目的，本发明采用如下技术方案。 

一种基于混合高斯模型的虹膜图像分割算法，包括内边界定位、外边界定位和眼睑检测； 

所述内边界定位具体流程如下： 

(S1.1)对输入的虹膜图像进行框架Laplace金字塔分解。 

(S1.2)基于混合高斯模型的瞳孔分割。 

(S1.3)基于形态学的反光区域消除。 

(S1.4)虹膜中心初始化。 

(S1.5)虹膜中心定位和内边界定位。 

作为本发明的进一步改进，所述虹膜中心定位为基于弦长均衡策略的虹膜中心搜索，具体步骤如下： 

第一步：虹膜中心的水平方向搜索； 

第二步：虹膜中心的垂直方向优化； 

第三步：计算内边界半径； 

所述外边界定位是基于区域平均灰度差的曲线演化方法实现外边界的定位，具体步骤如下： 

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

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

(S2.3)在图像中，基于区域平均灰度差的曲线演化实现外边界的定位， 

(S2.4)将上述迭代得到的外边界映射到原始图像I_r，即为虹膜外边界。 

所述(S2.3)中，基于区域平均灰度差的曲线演化实现外边界的定位，其迭代过程为： 

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₀作为外边界定位结果。 

所述眼睑检测对上下眼睑定位采用相同的方法，即首先利用水平方向的序统计滤波和极大值检测定位上下眼睑，再通过二次曲线进行拟合。 

所述眼睑定位过程如下： 

(S3.1)对内外边界定位后的虹膜图像进行水平方向的序统计滤波，得到图像I_e，消除 睫毛对眼睑的遮挡； 

(S3.2)在垂直方向上对I_e进行差分，得到二值图像I_ed。 

(S3.3)对二值图像I_eb进行二次曲线拟合，得到最终的眼睑定位结果。 

本发明的有益效果： 

(1)本算法将混合高斯模型引入虹膜预处理，实现了对虹膜内边界(即瞳孔边界)的准确定位。虹膜内边界定位的难点在于瞳孔的检测。瞳孔在虹膜图像中的灰度水平最低，一般采取阈值策略进行检测。本算法利用混合高斯模型所具有的对不同数据的自适应统计特性，对原始虹膜图像进行分割，获得灰度最低的瞳孔区域；然后在瞳孔内部初始化虹膜中心，利用弦长搜索策略实现内边界定位。 

(2)本算法在外边界定位过程中，同样在混合高斯模型分割的基础上，通过使局部对比度最大化实现了对外边界的准确定位。 

附图说明

图1本算法流程图。 

图2基于弦长均衡的瞳孔中心定位。 

图3外边界定位。 

图4基于区域平均灰度差的外边界迭代。 

图5对CASIA-IrisV3-Lamp数据库的处理结果。 

图6对CASIA-IrisV3-Interval数据库的处理结果。 

具体实施方式

下面结合附图与实施例对本发明作进一步说明。 

本方案提出的虹膜分割方法包括内边界定位(用S1.1-S1.5表示)、外边界定位(用S2.1-S2.4表示)和眼睑检测(用S3.1-S3.3表示)三部分。图1展示了完整的算法流程，步骤1-5为内边界定位过程，步骤6为外边界定位过程，步骤7-9为眼睑定位过程，各个过程具体为：步骤1.Laplace金字塔；步骤2.GMM分割；步骤3.阈值化及去反光；步骤4.基于弦长均衡的中心定位；步骤5.内边界定位；步骤6.外边界定位结果；步骤7.上眼睑检测；步骤8.下眼睑检测；步骤9：分割结果输出。 

内边界定位。 

内边界定位是整个分割过程的关键。我们利用混合高斯模型实现瞳孔和内边界定位。具体来讲，当对虹膜图像进行混合高斯建模后，由均值和方差分割，而瞳孔通常属于灰度最低 的区域。但是，由于睫毛或眼睑遮挡，分割出的瞳孔形状通常不是完整的圆形，我们又利用弦长均衡策略搜索瞳孔中心。内边界定位具体流程如下： 

(S1.1)对输入的虹膜图像I进行框架Laplace金字塔分解^[17]，其一层分解公式用矩阵形式表示为 

$\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(1\right)$

其中H和G为共轭镜像滤波器组。由(1)式得到概貌图像Ic和细节图像Id。此过程见图1中的步骤1。Laplace金字塔分解可以有效降低数据量，提高运行效率。 

(S1.2)基于混合高斯模型的瞳孔分割 

设定分类数K的值，我们将得到的概貌图像Ic通过混合高斯模型进行建模，其中灰度最低的部分被设定为瞳孔。由贝叶斯判据获得瞳孔区域的阈值Th，则 

$\mathrm{Ic}\_\mathrm{pupil}(i,j)=\left(\begin{array}{ccc}1,& \mathrm{if}& \mathrm{Ic}(i,j)\ge \mathrm{Th}\\ 0,& \mathrm{if}& \mathrm{Ic}(i,j)<\mathrm{Th}\end{array}\right)---\left(2\right)$

其中Ic(i,j)为(i,j)处的灰度值，Ic_pupil(i,j)表示阈值化的二值瞳孔图像像素值。 

(S1.3)基于形态学的反光区域消除 

由于反光，二值化后的瞳孔区域依然存在亮斑。我们利用形态学闭运算将其消除，即 

$\mathrm{Ic}\_\mathrm{pupil}=(\mathrm{Ic}\_\mathrm{pupil}\oplus b)\mathrm{\Theta b}---\left(3\right)$

其中Ic_pupul为二值化的瞳孔区域，形态学算子为b，形状为圆，和Θ分别表示膨胀和腐蚀算子。由此小尺度的反光点被上述操作消除而大尺度的瞳孔形状则得以保留。 

(S1.4)虹膜中心初始化 

为确保内边界的准确定位，虹膜中心必须初始化于真实的瞳孔区域内部。在步骤(S1.3)完成后，我们对二值化瞳孔图像进行水平和垂直投影，分别得到水平投影曲线P_h和垂直投影曲线P_v，并将它们各自极小值对应的位置初始化为虹膜中心(C_h,C_v)，记为点D。 

(S1.5)虹膜中心定位和内边界定位 

在虹膜中心初始化完成后，考虑到瞳孔的上部边界经常被睫毛或眼睑遮挡，我们由初始化点D分别沿0、π和3π/2三个方向引出射线，与瞳孔边界的交点分别为R,L和M。平移中心点的位置，当弦长|DL|,|DR|和|DM|三者相等时，可以认为当前点D即为真实的中心点。基于弦长均衡策略的虹膜中心搜索见图2，基于弦长均衡的瞳孔中心定位：步骤1：水平 搜索；步骤2：垂直搜索。具体步骤如下： 

第一步：虹膜中心的水平方向搜索（设定一个较小的正实数ε和正整数步长k）， 

第二步：虹膜中心的垂直方向优化 

第三步：计算内圆半径 

update DR; 

$\mathrm{radius}=\left(\frac{\left|\mathrm{DL}\right|+\left|\mathrm{DM}\right|+\left|\mathrm{DR}\right|}{3}\right);$

输出：虹膜中心D_c和内边界半径radius. 

外边界定位。 

与内边界的强烈灰度变化不同，虹膜外边界通常表现为变化缓慢的弱边缘，故采用一般的边缘检测算子效果并不理想。考虑到虹膜的区域灰度变化，本算法提出基于区域平均灰度差的曲线演化方法实现外边界的定位(见图3，外边界定位.，步骤1.选择扇形区域，步骤2.按照距离中心的远近进行序统计滤波，步骤3.外边界定位结果显示）。具体步骤如下： 

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

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

$\overline{I\_\mathrm{rp}}(i,j)=\left(\begin{array}{cc}\min (I\_\mathrm{rp}(i,j)),& (i,j)\in {I\_\mathrm{rp}}_{\mathrm{upper}}\\ \mathrm{med}(I\_\mathrm{rp}(i,j)),& (i,j)\in {I\_\mathrm{rp}}_{\mathrm{middle}}\\ \max (I\_\mathrm{rp}(i,j)),& (i,j)\in {I\_\mathrm{rp}}_{\mathrm{bottom}}\end{array}\right)---\left(4\right)$

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

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

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.4)将上述迭代得到的边界C₀映射到原始图像I_r，即为虹膜外边界。 

以上内边界定位算法要求边界内外区域要有足够的灰度差异，但是实际上很多虹膜图像并不满足这一要求，这会使迭代陷入局部极值。由此我们采用了序统计滤波^[18]，其步骤为：首先，对步骤(S2.2)所得的图像I_rp均分三部分，虹膜外边界位于中段；然后对上段(虹膜区域)进行极小值滤波，对中段(外边界所在区域)进行中值滤波以保持边缘，对下段进行极大值滤波。上述分段序统计滤波有效增强了虹膜图像的灰度对比。此外，分段序统计滤波为非线性操作，会在不同滤波区域的交界处产生明显的边界效应，这对轮廓迭代是非常不利的。我们将初始轮廓设定为图像中段内距离上边界(中段与上段边界)若干像素的直线，此直线向原理中心的方向演化，这样会克服边界效应带来的不利影响。 

眼睑检测。 

眼睑遮挡是造成虹膜图像质量下降的重要原因之一。本算法对上下眼睑定位采用相同的方法，即首先利用水平方向的序统计滤波和极大值检测定位上下眼睑，再通过二次曲线进行拟合。眼睑定位过程如下： 

(S3.1)对内外边界定位后的虹膜图像进行水平方向的序统计滤波，得到图像I_e，由此可以消除睫毛对眼睑的遮挡； 

(S3.2)在垂直方向上对I_e进行差分，得到I_ed，I_ed(i,j)为其像素值，将I_ed中每一列的最大值记为1，否则为0，即 

(S3.3)对二值图像I_eb进行二次曲线拟合，得到最终的眼睑定位结果。 

在近红外波段采集的仅包含眼部区域的虹膜图像，以通用格式(如bmp,jpg,tiff等)存储在计算机内，调用本算法进行处理，再将处理结果用识图软件显示即可。 

所述混合高斯模型数值实现方法： 

若虹膜图像灰度等级为N(一般为256)，令x_j表示第j个灰度等级，h(x_j)表示原始虹膜图像中第j个灰度等级的概率，prb_i(x)表示第i类的概率密度函数，K为分类数，prb_i(x_j)表示第i类中第j个灰度等级的概率，nloglik和loglik为初始化参数，则混合高斯模型的数值实现方式如下： 

第一步：概率密度函数的初始化 

第二步：概率密度函数的演化 

参考文献： 

[17]M Do,M Vetterli.Framing pyramids[J].IEEE Trans.Signal Processing,2003,51(9):2329-2342. 

[18]G R Arce,M Tian.Order statistic filter banks[J].IEEE Trans.Image Processing,1996,5(6):827-837. 

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。