基于区域Hough变换的眼睑退缩诊断

胡 敏，梁 挺，张光华，黄宏程

(1.重庆邮电大学 通信与信息工程学院，重庆 400065；2.太原学院 计算机科学与工程系，山西 太原 030000)

0 引 言

甲状腺相关性眼病(thyroid associated ophthalmopathy，TAO)是一种与自身免疫有关的眼部疾病，是一种常见的眼部疾病，占Graves病的40%～50%，免疫性甲亢的2～5%。关于眼睑退缩等定义见参考文献[1-3]。目前已知的眼睑与虹膜外圆的位置信息可分为全局和局部两种。

目前国内外有关虹膜位置的相关文献多以某一特定的虹膜采集设备为基础，收集到的图像如下：①仅包括人类的眼睛；②对虹膜及瞳仁进行了扩大；③增加了对眼、眉等的干扰作用；④虹膜、瞳孔和巩膜有显著的区别。传统的、通用的人脸图片都是由便携式设备获取的，例如手机和普通相机，而不是由仪器获取，已有的研究手段无法对其进行有效应用，因此，开展基于面部的虹膜定位的技术十分重要。人类的虹膜和瞳孔都位于眼眶的正中央，上眼睑的形状和位置都有相同的特点，而瞳孔、虹膜、眼睑等都比正常的瞳孔、虹膜、眼睑等情况更加复杂。首先，要解决TAO病人影像上的虹膜位置问题，然后将此算法应用到便携式上。

全采集采用便携式设备，通常为全貌、局部采集使用专业的虹膜仪器，仅含有眼睛的部分信息。虹膜的定位一般是确定虹膜的内外边缘和上、下4个分界线[4-6]。基于全局采集的虹膜定位算法包括：Daugman等[5]采用积分微分算子，通过圆型模板对虹膜的中心及半径进行搜索，得到了较好的定位效果。Wilde等[7]提出了利用霍夫变换进行虹膜内、外边界的方法，但是要对整个图像进行检索，不仅费时费力，而且耗费较大。目前，基于局部提取的虹膜定位方法主要有：田子林等[8]所提出的基于最小二乘法的虹膜定位方法；王琪等[9,10]所提出的基于SDM和SIFT的虹膜定位算法，通过对上睑和虹膜外缘的比例变换不变特征进行定位，从而使其定位精度得到了提高，但其计算量大、耗时长。采用特殊的摄影器材进行局部采集，可以放大眼睛的特点，但由于设备的成本较高，不宜大规模购买；而且，近距离获得的眼球信息不能达到眼睑退缩的判别准则。

为了解决上述问题，本文在利用便携式摄像机采集到的完整的面部影像基础上，采用区域霍夫变换进行虹膜定位，并与EUGOGO眼睑退缩准则相结合。实验结果表明，该算法能够很好地实现对图像的辅助诊断。

1 区域Hough变换的虹膜定位方法

1.1 特征点提取

样本图像由便携式设备拍摄，以第一眼位图为标准(当两眼在真正水平面注视6 m以上的目标时，两眼注视线同时向前方相并平行)，背景颜色选择白色，图像包含了从头发到双肩之间的整个面部区域。面部对齐是根据图像中的区域和尺寸来确定面部形状[11]。脸部特征点是面部对齐的重要指标，被看作是学习一个回归函数的过程[12,13]，将原始图像作为输入，利用多个简化的小波变换算法对所抽取的特征进行分类，目标的输入取决于上一步的输出，反复地进行残差拟合直至获得最优的估计形态。初始的特征点定位共包含68个特征点，并没有覆盖额头区域，如图1所示。

图1 68特征点效果

在此基础上，本文利用拓展后的81特征点法，如图2所示，该方法在前额头上新增13个特征点，弥补了68个特征点的不足之处。所有的特征点依次由1-81标记，图2中1至17的标号与69至81的标号，30个特征点共同组成面部外轮廓，眼睑部位一次标号为37至48。

图2 81特征点效果

1.2 Bezier曲线拟合

本文采用Bezier曲线拟合面部轮廓特征点，截取面部图像，排除头发、背景等干扰因素。Bezier曲线是法国工程师Bezier[14]在1962年为了设计汽车车身形状提出的，应用于二维图形程序，在计算机图形学中是相当重要的参数曲线。与最小二乘法等拟合方法相比，其优势在于Bezier曲线实现简便、需要的条件更加简单且拟合后的曲线更光滑。

贝塞尔曲线的定义严格依赖于确定该段曲线控制点的个数，控制点也就是曲线上的点，N次多项式的曲线由N+1个顶点确定。贝塞尔曲线上各点的参数方程为

(1)

其中，Pi为第i个顶点的坐标值，基函数Bi,N(t) 为Bernstein多项式，该多项式公式为

(2)

由上面两式可得Bezier曲线，根据N+1个顶点 (A0,A1,…,An)， 得到其贝塞尔曲线公式为

(3)

该公式的原始数学表达式为伯恩斯坦多项式，即在 [A0,An] 区间上所有的连续函数都可以用多项式来逼近，并属于一致收敛。曲线的起点是A0，终点是An,在所有控制点都是同线的情况，曲线表示为一条从A0到An的直线。该曲线可以在任何控制点上进行分割，每条被分割的曲线都属于Bezier曲线，其两端的点分别为分割处的控制点。

在1.1节中，图2通过标号1-81标注了所有的特征点，面部最外一圈的特征点构成了拟合区域，包含了标号1-17，69-81，一共30个特征点。将30个特征点用A0-A29表示，即用于人脸特征点拟合的Bezier公式为

B(t)=P0(1-t)30t0+30P1(1-t)29t+…+
30P29(1-t)t29,t∈[0,1]

(4)

通过式(4)的拟合效果如图3所示，图3闭合线条区域为30个特征点拟合后的效果，该区域切割后的效果如图4所示。对比图3和图4可知经过Bezier曲线拟合后头发、双肩以及左右耳朵等明显的干扰区域被处理。

图3 Bezier曲线拟合效果

图4 处理后效果

1.3 三次B样条插值

本文应用B样条插值方法[15]获得了眼睑轮廓曲线，并充分利用了它的局部支撑性和插值足够光滑的特性优势，从而获取了更加精准的上下眼睑轮廓曲线。B样条曲线包含一系列的控制点，并将其划分为节点和系数，每一个系数对应一个控制点，其核心思想是采用分段低阶多项式将高次多项式进行连续拼接。由于Bezier曲线不能进行局部修改，而且在拼接成长线时，会产生一些波形，对结果产生影响，因此，我们选择了B样条进行轮廓拟合；而上下眼皮分别用4个点来进行定位，它们符合三次B样条曲线的对应条件，也就是三次B样条曲线的插值。

在对上下眼睑进行划分时，左右两个眼部被均分为上下两个部分，包括4条眼睑线，分别称之为左上眼睑、左下眼睑、右上眼睑、右下眼睑。在1.1节图2中，可将4条眼睑分别对应4个标号，即三次B样条插值所需的4个控制点。左上眼睑对应的4个特征点标号为：(37、38、39、40)。左下眼睑对应的4个特征点标号为：(37、42、41、40)。右上和右下分别对应的标号为：(43、44、45、46),(43、48、47、46)。

B样条曲线的总方程为

(5)

其中，Pi是控制曲线的特征点，Fi,k(t) 则是K阶B样条基函数。

三次B样条曲线方程中基函数为

(6)

(7)

(8)

(9)

(10)

由式(6)到式(10)可得三次B样条曲线方程

P(t)=P0*F0,3(t)+P1*F1,3(t)+
P2*F2,3(t)+P3*F3,3(t)

(11)

下面是三次B样条曲线的拟合法，如图5所示，图5中封闭线条区域表示上下眼睑经过式(11)插值后的效果区域。结果表明，采用三次B样条插值算法进行上下眼睑边缘定位有很好的效果，同时有很高的准确率。

图5 眼睑拟合

1.4 Hough圆检测

Hough变换是Paul Hough在1962年首次提出的一种线条检查算法，在1972年，经过Richard Duda和Peter Hart推广使用，经典的Hough变换通常被用来检测图像中的直线，经过众多学者的不断研究和发展，Hough变换被扩展到可以进行任意形状物体的检测，现应用最为广泛的形状是圆和椭圆的检测[16,17]。人眼的虹膜区域通常位于巩膜和瞳孔之间，其外边界和巩膜相切，实际为眼白和棕色线。虹膜在形状上无限接近于圆形，并且具有不会轻易受到其它因素影响而改变形状的特性。针对虹膜这种特性，本文利用Hough变换圆检测方法进行虹膜外圆的定位。

Hough变换圆检测原理是对于一个半径为r，圆心坐标位置为 (a,b) 的圆，将其表示为

(x-a)2+(y-b)2=r2

(12)

圆检测方法步骤如下：

(1)建立一个由各个象素单位组成的累积空间。开始时，将每一个单元格都设定为0。

(2)针对各图象的边界点座标位置(x,y)，可以是一个圆形的中央的网孔数值，如方程式(12)所示。方程式中的单位用a代表。

(3)在上述的步骤中，从各个可能发现的数值a中，得到符合方程的全部可能的数值b。

(4)在累积空间内搜寻局部极大值。这些单元格表示算法检测到的所有圆圈。

以上所述的圆检法所涉及到的参数空间是三维的，若要在三维空间上进行证据累计，则需要耗费大量的时间和空间复杂性。因此本文对该Hough圆检测方法进行改进，主要考虑优化空间开销和时间开销。为此，对这种Hough圆法进行了改进，重点是在优化的空间和时间上。本文提出了一种基于局部霍夫圆法的新算法，该算法可分为两个阶段：

(1)选择目标区域。基于Bezier曲线方法，拟合图像，然后根据人脸五官的分布，获取目标眼部区域信息。

(2)确定搜索区域范围。首先，利用改进后的Canny算子对图像进行边缘信息提取，然后对搜索区域进行去噪处理，这样更好地突出虹膜外圆的线条信息，并规定了经典Hough圆检测方法中的搜索半径r指定范围值，而非盲目进行区域搜索，避免了三维的参数空间，只需要在目标区域内对检测圆心(a,b)进行证据累计，从而缩短了时间和空间，增加了虹膜外圆定位的准确率。

目标区域选择通过缩小Hough圆圈来探测和降低空间消耗。该算法采用了一种改进的Canny边界提取算法来突出图像的边缘，同时根据虹膜的直径来限定搜索半径r的范围，从而降低了搜索的耗时。

1.4.1 选择目标区域

“三庭五眼”的概念最早在明末王绎的《写像秘诀》中记录[18]。“大三庭”通常是把头部上下分为三等分，从前额发际线至眉骨，从眉骨至鼻底，从鼻底至下颏，各占脸长的三分之一。“五眼”从左侧发际至右侧发际，为5只眼形。两只眼睛之间有一只眼睛的间距，两眼外侧至侧发际各为一只眼睛的间距，各占比例的1/5。如图6所示。

图6 “三庭五眼”

本文截取眼部区域位于图6中的“中庭”处。直接取“中庭”区域为眼部区域会包含部分眉毛信息，为了避免其干扰，本文通过大量的图片测量和验证，确定目标区域为眉线以下六分之一以及鼻底线以上二分之一。同时对目标区域进行归一化处理，统一像素值为278×59。以4位TAO患者图像为例，截取目标区域图，分别如图7(a)、图7(b)、图7(c)、图7(d)所示。

图7 TAO患者目标区域

1.4.2 确定搜索区域范围

搜索区域处理主要分为两个部分：

(1)区域圈定：Canny算子能检测真正的弱边缘，可使用两种不同的阈值分别检测强边缘和弱边缘，不容易受噪声的干扰，具有准确的边缘定位功能，所以本文采用Canny算子对得到的目标区域进行边缘检测。经过Canny算子处理后的图像存在无用的小面积区域，本文将结合连通域处理方法进行图像筛选，过滤无用区域。改进后的Canny算子工作流程如图8所示。

图8 改进后Canny算子工作流程

图9 原始Canny算子处理后效果

图10 改进后Canny算子处理效果

图9(a)、图9(b)、图9(c)、图9(d)这4幅图像经过原始Canny算子处理后眼睑与虹膜基本轮廓较为明显，每位患者的左右眼虹膜区域内存在少数小圆圈(小圆圈是拍摄时人眼部位的反光聚焦点，表现为白色亮斑，经Canny算子处理后表现为小圆圈)，眼睑处存在断断续续的不连接线段。由于TAO患者眼部区域的特殊性和每个人眼部区域的差异性，部分图像经过Canny算子检测后虹膜轮廓不能呈现出完整的轮廓信息。图10(a)、图10(b)、图10(c)、图10(d)是经过改进后Canny算子处理结果图。对比图9和图10，小圆圈和不连接线段均被去除，保留了大部分虹膜外圆的边缘信息。

(2)区域搜索：由于人类的虹膜直径比较相近，平均虹膜直径12 mm，而中国人的平均虹膜直径大约为11.4 mm。本文对随机的500张TAO患者图像和100张非TAO患者图像进行虹膜直径的估计，98.6%的TAO患者眼直径在10 mm～13 mm范围内，1.4%的TAO患者眼直径在该范围外，所有的非TAO患者图像眼径在10.5 mm到11.6 mm之间。直径分布见表1。

表1 直径分布

根据直径分布表信息，1000张TAO患者图像中只有14张图像的眼直径范围在10 mm～13 mm外，即将搜索半径r的值缩小在10 mm～13 mm内。同样以4位TAO患者图像为例，原始Hough圆检测实验效果如图11所示，图中4位患者的定位结果存在两个主要问题，第一，检测到左右眼虹膜外圆但外圆的大小和位置相比于真实的虹膜外圆位置和大小都出现了一定偏差；第二，部分图像只能检测左眼或者右眼的虹膜外圆轮廓。本文的区域Hough圆检测实验效果如图12所示，图中本文方法对每张图都能准确找到两只眼的虹膜外圆，保证虹膜外圆大小和位置符合其实际值，另外对虹膜外圆边缘信息不足的情况，补全了拟合圆的轮廓信息，完成虹膜外圆的定位。

图11 经典Hough圆检测定位结果

图12 本文Hough圆检测结果

2 实验结果和分析

本文实验的数据集由两部分组成，分为私有数据集和公有数据集。私有数据集主要用来完成诊断眼睑退缩实验，由眼科医院提供的2500张图像组成，其包括1822张已知存在眼睑退缩的TAO患者图与678张确认为非TAO患者的图像。图片尺寸统一为1915×1080，图片为jpg格式。公有数据集主要用来验证本文所提出的方法同样适用于专业设备拍摄的局部图像，在定位精度和时间上有一定的提升，数据集从中科院CASIA虹膜库中随机抽取的100幅图像组成，分辨率为640×480，图片为png格式。

2.1 精度分析

为了进一步验证所提出的方法的可行性，进行了两组对比实验，实验基于中科院CASIA虹膜库随机抽取的2500幅图像。第一组是通过对虹膜的定位精度、虹膜外圈的定位时间和最小二乘法与霍夫变换的虹膜定位算法相比较，后者采用形态学、自适应阈值算法、Canny边界提取、最小二乘法等技术，对虹膜外圆进行了定位。第二组是将本文方法与文献[9]所提出方法进行对比，后者根据SIFT技术和SDM技术，结合文献[9]中所提的虹膜定位技术，采用径向对称变换、微积分算子、SIFT特性以及SDM等方法，实现了对目标的精确定位。并给出了式(13)中的位置精度的计算公式

(13)

式中：O表示实际虹膜面积，N表示测得的虹膜面积，C表示重复的虹膜面积。根据实际虹膜与测得的虹膜面积覆盖率表示其定位精度。实验结果见表2。

表2 CASIA虹膜库图像定位结果对比

表2在CASIA虹膜库数据集的基础上一共进行了5组实验，根据欧洲眼病专家组(European group of Graves’ orbitopathy，EUGOGO)[19]进行精度和方差的对比，Daugman提出的基于微积分的虹膜定位算法、Wilde提的Hough变换与Canny相结合的定位算法、文献[8]中的最小二乘法与Hough变换相结合的虹膜定位算法、文献[10]提出的基于SDM和SIFT虹膜定位算法以及本文提出的基于区域Hough变换的虹膜定位方法。实验数据集从CASIA数据库中随机抽取2500张，从表格可以看出，定位精度、方差和外圆定位耗时三方面进行对比，本文所提方法在精度、稳定性以及时间上都优于其它算法。其中Daugman和Wilde所提的定位方法由于在早期就已经提出，所以其精度较其它算法都较低。文献[8]中二值化过程比较粗糙，不够精确，其利用灰度直方图的波峰波谷取得阈值，很容易导致实验结果出现一定偏差，但是在区域Hough变换法中就通过对目标区域和搜索半径的范围进行确定，在基于CASIA虹膜库数据集的实验中表明其具有更高精度和更低耗时的优势。文献[10]中方法主要思想是以径向对称和微积分算子为基础确定虹膜中心和半径，但是对图像质量要求较高，初步的计算结果将直接影响最终的虹膜定位结果。综上所述，本文所提区域Hough变换法进行虹膜定位在较短的时间消耗情况下具有较高的精度。

2.2 眼睑退缩诊断结果分析

为了验证本文方法在实际应用场景中能有效辅助医生诊断患者是否存在眼睑退缩，共进行了5组对比实验，同时计算了实际误诊率。在实验中对2500张图像进行眼睑退缩诊断，其中1822张确认为TAO患者，678张确认为非TAO患者图像。最终实验定位效果如图13所示(图片已经获得可公开授权)。

图13 虹膜外圆与眼睑定位

眼睑退缩诊断结果见表3。

表3 眼睑退缩检测实验结果

如表3所示，进行了5组实验，在2500张数据集上通过准确率、误诊率和平均耗时3个指标来评估算法。其中Daugman和Wilde这两种方法相比其它算法存在准确率低且误诊率高的问题，这是因为这两种算法对图像质量要求较高，但是患者图像面部的特征不清晰并且不规则，所以在检测时存在较大误差，耗时长。文献[8]和文献[10]准确率与误差值较为相近，耗时也相近，这是由于两种算法均是在原始方法上进行部分特征改进，所以在准确率和误诊率以及耗时上比Daugman和Wilde方法更好，但是由于患者图像的特殊性，其准确率等都存在一定上限，而本文的算法是专门针对固定的患者图像提出的，所以在准确率等指标上比其它方法更具有优势与针对性。本文算法准确率为93.5%，且误诊率低，达到了实际应用的基本要求。综上所述，本文所提出的区域Hough法在准确率、误诊率以及耗时三方面均优于其它几种算法。

3 结束语

针对传统的人工眼睑退缩诊断准确率低以及耗费时间长的缺点，本文提出了一种基于区域Hough变换新的虹膜定位方法。与传统的Hough方法相比，该算法采用特征点、Bezier曲线和三次B样条拟合出眼睑轮廓和目标检测区域，从而缩小了搜索范围，并减少大量的时间开销，通过改进的Canny边缘检测方法和改进的Hough圆检测方法提升了巩膜外圆定位的精度，再结合EUGOGO眼睑退缩标准进行诊断。实验结果表明本文改进的方法对眼睑退缩的诊断有较大的提升，满足实际计算机辅助诊断的需求，又能在一定程度上有效提高在利用便携式设备局部获取图像下的精度，缩短定位时间[22]。