生理结构先验引导下的虹膜精确分割算法

年炳坤，丁建睿，史梦蝶，黄子晨

(哈尔滨工业大学(威海)计算机科学与技术学院,山东 威海 264209)

随着生物识别技术的发展，人们逐渐将某些具有区分能力的生物特征运用到身份识别中。其中虹膜是所有生物识别系统中应用最为广泛的生物特征。虹膜是位于黑色瞳孔和白色巩膜之间的圆环状部分[1]，拥有高于人脸和指纹识别几个数量级的精确性，被广泛应用在身份识别的各个领域，如：美国-墨西哥边境控制，2004年雅典奥运会身份鉴定，阿联酋驱逐人员的管理等[2]。

20世纪90年代Daugman[3]在假设瞳孔与虹膜为同心圆的前提下提出了虹膜分割算法，取得了一定的效果。而后针对各种不同的情况，许多学者针对某一种特殊情况(瞳孔的中心偏离[4-5]、遮挡[6-7]、镜面反射[6])，提出了特定的分割方法。为了进一步提高分割准确度，Daugman[8]又提出了基于傅里叶变换并假设瞳孔与虹膜为非同心的Snake主动轮廓模型，随后Sutra等[9]彻底抛弃掉了同心以及轮廓外形为圆的假设，提出基于Viterbi的主动轮廓探测算法，可以较为精确地探测虹膜与瞳孔，但是算法对非虹膜区域(睫毛、眼睑)的干扰较为敏感。Llano等[10]利用主成分分析(principal component analysis,PCA)进行图像融合，随后提出了基于拉普拉斯金字塔的融合方式(modified Laplacian pyramid,MLP)[11]，进一步提高了虹膜系统的性能。但Daugman的算法将瞳孔与虹膜看作同心圆，这就会导致损失大量的有效信息；而Snake主动轮廓模型和Viterbi算法虽然抛弃了非圆、非同心的假设，但是算法本身受睫毛、眼睑等噪声影响较大；而PCA融合算法和MLP融合算法虽然有效地提高了系统的识别性能，但其本质上是对已有分割算法进行融合，让分割结果更容易被特征提取和编码，因此其性能受分割算法的制约。

同时国内学者刘琳等[12]对于主动轮廓算法进一步研究，提出了结合KL散度和RSF模型的主动轮廓图像分割方法[12]，提高了图像分割的精度，可以一定程度上提高识别精度。

本文基于CLBP(completed local binary patterns)思想[13]，提出一种形状敏感的成分分析算法，可以有效检测虹膜区域，同时减少其他成分(如：睫毛与眼睑区域)的干扰，从而避免分割过程中虹膜的信息丢失，提高整体的分割精度，进而提高虹膜识别系统的性能。

1 本文方法

本文提出的算法由粗分割与精细分割两部分组成，整个流程见图1，在分割基础上利用Daugman[13]提出的方法进行匹配。

图1 本文方法流程图Fig.1 Flow chart of proposed method

首先利用瞳孔灰度值低于其他区域的生理特性，对图片进行二值化。直接使用二值图像进行后续处理存在以下问题：1)二值化图像存在噪声点，并不是所有像素点都属于瞳孔，见图2(a)；2)由于拍摄条件和瞳孔本身具有反光的生理特性，属于瞳孔的像素点在局部区域会出现异常，致使二值化图像中瞳孔可能出现断裂现象，不利于后续处理，见图2(b)。

由于噪声点大多孤立，或者小部分局部聚集；同时瞳孔局部灰度值异常面积不会太大。为解决上述噪声问题，本文对二值图像进一步进行形态学操作，去除噪声的同时，可以在一定程度上对瞳孔断裂区域进行连接。形态学腐蚀、膨胀、闭操作分别见式(1)～(3)[14]：

A⊖B={z|(B)z∩Ac=∅}

(1)

(2)

A·B=(A⊕B)⊖B

(3)

式中A为输入的二值化图像，B为结构元。由于本文所处理的形状多为圆形，所以选择3*3的结构元。处理结果见图2(c)。

图2 二值化图像示例Fig.2 Examples of binary images

1.1 瞳孔粗分割

通过上述操作，可以认定图像中的像素值为0的点绝大多数属于瞳孔。先假定瞳孔为圆，首先利用式(4)计算图像中黑色点集的几何中心，如下所示：

(4)

式中：I(xk,yk)为坐标点(xk,yk)的像素灰度值，N为图像中灰度值为0的像素个数，(x,y)为黑色点集的几何中心，为瞳孔内一点。

对具有最大面积的连通区域补全后，求得其最小外接矩形，该矩形的内接圆边界即为瞳孔的粗边界，该边界用于补全晶状体、眼镜等反光介质导致的局部灰度值异常所带来的瞳孔断裂问题。

1.2 瞳孔精确分割

利用连通域分析，找到最大连通域，并标记出其边界，这是瞳孔的精确边界，记作A(I)。进一步采取如下规则将瞳孔粗边界R(I)和瞳孔细边界A(I)进行融合，以解决由于瞳孔自身局部灰度值异常可能导致的问题。

1)如果粗边界完全包含在精确边界之内，则采用精确边界；

2)如果粗边界部分包含在精确边界之内，则最终边界为粗边界和精确边界的并集。

其数学描述如式(5)所示：

F(I)=

(5)

式中：I代表输入的原始虹膜图像，R(I)为瞳孔粗边界，A(I)为瞳孔精边界，F(I)代表融合后的掩码图，融合后的边界即为瞳孔的最终外边界,结果见图3(a)。

1.3 虹膜外圆粗分割

在假设虹膜外边界为圆形的分割算法中，大多采用Daugman的球形积分算子，在理想情况下，即没有眼睑、睫毛遮挡，虹膜完整呈现在眼球当中时，该方法的准确度较高，如图3(b)-(1)所示。

但是大多数情况下，虹膜图像大多存在遮挡问题，此时，球形积分算子达不到预期效果，如图3(b)-(2)所示。

经过分析，影响精度最大的因素是上眼睑的存在。另外，球形积分算子本质上是探测梯度变化，当检测到眼皮和眼球分界线时，边界线灰度值变化远比虹膜和巩膜分界线的要大，所以影响较为显著，从而导致算子检测错误，如图3(b)-(3)所示。

图3 过程化分割结果示例Fig.3 Examples of procedural segmentation

针对以上问题，本文对其进行改进。根据眼球的生理特性可知：眼睑、睫毛的遮挡只会发生在瞳孔的上下区域，而左右区域一般不会被遮挡，只存在巩膜和虹膜的分界，而这有利于我们的探测，利用此先验知识为Daugman的球形积分算子增加约束，得到新的弧形积分算子，如式(6)所示：

(6)

式中：r为所求虹膜的半径，(x0,y0)为瞳孔几何中心的坐标，I(xk,yk)为图像中坐标为(xk,yk)的像素点的灰度值。

利用式(6)可以得到粗分割的虹膜外边界,同时可以有效解决不同程度的遮挡问题，结果见图3(c)。

1.4 虹膜外圆精分割

根据CLBP中的思想[13]，一幅图像可以初步处理成dp，将原图像划分成多个3*3的区域集合,dp的定义如式(7)所示：

dp=gp-gc

(7)

其中gc为3*3区域的中心点，gp为其相邻的点，而dp可以进一步被分解为两部分，sp和mp，定义如式(8)所示：

(8)

文献[13]中同时证明了mp中含有图像绝大部分信息，而sp中含有较少的信息。本文认为：mp中的信息为结构信息，而sp中的信息为关系信息。为充分利用两种信息，提出以下规则对其进行融合，起到类似滤波的效果，但比传统意义上的滤波更加具有导向性，能够在探测指定形状的同时，减小其他噪声成分的干扰，同时可以提高精度。

式(1)中,Ai为不同浓度样品与DPPH溶液的混合液的吸光度；Aj为不同浓度样品与无水乙醇的混合液的吸光度；A0为DPPH溶液与无水乙醇的混合液的吸光度。

1)如果需要检测的形状为横向长条形，则采用上下4连通域，如果为竖向长方形，则采用左右2连通域；

2)如果需要检测的形状为圆形，则采用8连通域；

3)如果想扩大搜索区域，则可以扩大连通域。

1.4.1 睫毛分割

睫毛不同于眼睑，可以近似看成一个各向异性的长方体，由此可以构造一个对四周敏感的融合式。本文构造了一个8连通域区域，中心区域的坐标点为(0,0)，集合G={(1,0),(0,1),(-1,0),(0,-1),(1,1),(1,-1),(-1,1)，(-1,-1)}。为了减少计算量，本文定义先前虹膜外圆粗分割时得到区域为ROI区域，在此区域内进行运算。对ROI区域中所有c∈C和ε∈G使用式(9)、(10)

(9)

(10)

得到滤波后的结果作为睫毛的掩码图层，见图3(d)。

图3(d)-(1)证明本文提出的算法可以检测出较为稀疏的睫毛，而图3(d)-(2)和图3(d)-(3)则证明算法对于密集的睫毛同样具有很好的检测效果。在ROI区域内，无论是上眼睑的睫毛，还是下眼睑的睫毛，都被完整的检测。同时，结果也符合先前的构造规则。

采用较高运算量的长条状8连通域进行检测，检测结果基本都是长条状物体，而上下眼睑边缘处阴影由于拍摄原因，可以近似认为是黑色的不连续阴影，算法也可以准确检测，这进一步证明了式(9)和式(10)对于给定形状的敏感性和稳定性，能够较大程度上去掉其他噪声点干扰。

通过虹膜外圆粗分割，可以得到虹膜ROI区域，但是这个ROI区域包含了大量不需要的信息，例如：眼球以外的眼睑区域。探测出整个眼皮区域的难度很大，可以探测眼球和眼睑的边界区域，然后和虹膜和巩膜边界一起构成一个闭合区域，这个区域就是去除眼睑后的掩码图层。

从灰度值来看，这个边界区域不一定是连续的，可能存在眼泪等反光物质，在边界线的某一点处，存在灰度值异常，从而断裂开来；同时眼睑的探测还需解决如何在睫毛的干扰下，有效探测出准确边缘。

根据眼睑的生理学特点可知，人的眼睑边界可以近似看成一个横向长条状，同时为了和睫毛区分，提高探测精度，本文选择上下4连通域，只探测横向的相关区域，来屏蔽竖向睫毛的干扰。构造一个4连通域，连通域中心点的坐标为(0,0)，G={(0,1),(0,2),(0,-1),(0,-2)}。对所有c∈C在ROI区域中的像素点使用式(9)和(10)。

通过连通域分析后，结果见图3(e)。可以看出，新构造的算式，能够准确地探测出眼睑的边界区域，如图3(e)-(1)结果所示，边界上存在局部反光导致的灰度值异常区域可以被准确避开，最终形成一个闭合区域；图3(e)-(2)中，在如此多睫毛的干扰下，仍能准确地检测出边界；图3(e)-(3)中，算法可以有效避开睫毛密集的地方，虽然损失了睫毛间隔之间小部分属于虹膜的区域，但成功避免了大量睫毛带来的误差。

最后将多个掩码图进行合并，可以得到最终的虹膜分割结果，对其进行特征编码并识别。

2 实验设计

本实验采用开放数据集CASIA-IrisV3-Interval和CASIA-IrisV1[15-16]。虹膜区域分割后，采用OSIRIS 4.1版[17]的匹配模块进行识别，它是在Biecure项目框架下的一个开源虹膜识别系统。

两种数据集中照片格式均为BMP，同时数据集中的图像质量并不相同。本文所提出的算法采用的是传统分割技术，不涉及到机器学习内容，不存在过拟合等现象。

CASIA-IrisV3-Interval数据集中包含了来自247个不同的个体，共计2 639张图片；CASIA-IrisV1数据集中包含了来自104个个体，共计709张图片。数据集中的图像示例见图4。

图4 虹膜数据集示例Fig.4 Examples of iris datasets

考虑到虹膜识别本质上是一个分类问题，本文采取衡量分类性能的指标，PGAR(正确接受率)在PFAR(错误接受率)为10-4时的值,PEER(等误差率)以及PFRR(拒识率)的值，如公式(11)～(13)所示。对于一个二分类模型的单个样本一共有4种预测结果：RTP(true positive)、RFP(false positive)、RFN(false negative)、RTN(true negative)。

(11)

(12)

(13)

PEER值为PFRR与PFAR相等时的值。同时采用ROC曲线下的面积(AAUC)来衡量算法的分类性能。

3 实验结果

图5选取了对比实验中涉及到分割的算法：CHT[18](以同心圆为假设并部分改进的算法)、Viterbi以及本文所提出算法的内外虹膜分割结果样例。然后将分割结果基于Daugman提出的匹配方式输入到OSIRIS Version 4.1[15]系统中，得到任意两幅图片的相似程度，然后设定阈值，计算出本阈值下的PGAR、PEER值，并求出对应的AAUC值，实验结果见表1。

从表1中可以看出，CHT分割算法的性能最低，原因是这个算法基于圆假设，包含了一部分无用区域或者少包含了一部分有用区域，如图5(b)-(2)，去掉了眼睑遮挡区域的同时，也去掉了很大一部分有用的虹膜区域，如图5(c)-(2)中所示。Viterbi算法很大程度上提高了准确性,对于数据集CASIA-IrisV3-Interval，Viterbi算法将PGAR在PFAR=10-4的值从95.32%提升到了97.04%。这表明Viterbi算法具有更好的分割性能，提升的原因在于Viterbi算法采用了主动探测轮廓，能够较为精确地定位出瞳孔、虹膜的外边界区域，而不是宽泛的规定和去掉某些区域，从而有效区域损失和无用区域的包含远比CHT算法要小的多，但Viterbi算法也存在缺点，如图5(d)-(3)中显示，算法对于眼睑边界部分检测不够准确，当边界不太明显且存在局部灰度值异常的时候，存在较大范围误差，同时如图5(c)-(3)图像中所示，由于假设分割区域灰度值成正态分布，在给定2.35个σ(标准差)外都认为是杂质，并加以去除，所以导致丢失了部分虹膜信息，影响了分类性能。

表1 在数据集CASIA-IrisV3-Interval上的结果Tab.1 Results on database CASIA-IrisV3-Interval

图5 分割结果示例Fig.5 Examples of segmentation results

基于PCA主成分分析和MLP的融合算法，融合了CHT和Viterbi算法的分割结果，一定程度上弥补了缺陷，但是受CHT影响较大，并不能很好地提高算法性能，甚至有可能会下降。PCA主成分分析的原理是利用数学上的PCA降维手段，将两个对象的数据，转换到K个特征向量构成的新空间中，而MLP融合算法采用拉普拉斯金字塔进行融合，本质上合成新的特征向量，这就导致融合结果强烈依赖于融合对象，也就是分割算法的性能，如果两个分割算法能够互相弥补自身不足，那么融合表现就好，反之如果缺陷都较大，那么融合算法一定程度上也会放大缺陷。如在表2中，PCA融合后的PEER值反而上升。而本文提出的算法对于同种结构的图片，有较强的鲁棒性，可以根据不同目的构造不同的算子，将睫毛、眼睑等虹膜特有噪声检测出来，在精确分割的同时，减少虹膜信息的损失。

表2 在数据集CASIA-IrisV1上的结果Tab.2 Results on database CASIA-IrisV1

从表2中可以看出，表现最好的MLP融合算法和本文的算法结果相差不大，PGAR、PEER值基本一致，而AAUC值也只是略高一点。通过分析CASIA-IrisV1数据集的特点，得出原因在于数据集CASIA-IrisV1拍摄条件较好，数据集中的虹膜图像，无论是眼睑、睫毛等遮挡均比CASIA-IrisV3-Interval中要少，图片整体质量较高。所以Viterbi和CHT算法的缺陷(对于睫毛、眼睑的探测精度不足)被弥补。PCA和MLP融合算法融合了Viterbi和CHT的优点，一定程度上特征被增强，使得融合后的图像更利于分辨。导致本文提出的算法在新性能上并没有显著的提高。同时这也说明了本文提出的算法，在分割性能上具有优越性，在没有融合的情况下达到了相似的性能，这也意味着如果将PCA与MLP融合算法的融合对象CHT换成本文算法，那么精度会进一步提高。

算法在保证精度的同时，时间复杂度较低。算法中的大部分运算可以表示为简单的加减运算。这大大提高了算法的性能。表3给出了本文算法与CHT和Viterbi算法的运行时间比较。这些算法在带有Intel(R)Core(TM)i7-7700hq、16 GB内存、GTX1060显卡的计算机上运行。在Windows 10上使用C++编程语言和OpenCV4.1.0。

表3 算法平均运行时间Tab.3 Average running time of different algorithms s

4 结 论

本文基于CLBP思想，提出了一种成分分析算法，并且在开源数据集CASIA-IrisV3-Interval与CASIA-IrisV1上进行了实验，与CHT、Viterbi、PCA和MLP融合算法进行了比较，从实验数据上看，本文所提出的算法，在两个数据集上都具有较高的AAUC值和较低的PEER值,说明本文算法具有更优秀的分割性能。

同时，在实验过程中发现本文所提算法对于低质量的图片，如有模糊、散焦等问题的图片分割效果略差，并且算法较强依赖于进行二值化的阈值，这导致算法对于不同摄像设备所拍摄的不同结构的虹膜图像效果不显著。如何根据图像质量自动化确定合理的阈值，是今后算法改进的主要方向，可以借助深度学习方法，自主决定阈值，来提高算法的鲁棒性。