一种基于图像处理的人眼全自动角膜曲率计研究

赵俊奇 郭智勇 陈安世 刘海峰

1(中北大学光电仪器厂，太原 030051)

2(清华大学环境工程系，北京 100084)

引言

角膜是眼球前面一层透明组织，是眼屈光系统中最大的折射面，有一定的曲率半径和屈光率，人眼的全部屈光率是58.64 m－1。在眼的屈光度范围中，很大部分的作用是由角膜表面完成的。角膜表面屈光率是40.0～45.0 m－1，人眼角膜屈光参数包括曲率半径和角膜屈光度(包括角膜球镜、柱镜和轴角)。隐形眼镜的科学验配、眼睛角膜屈光度的矫正手术、对某些角膜病(如圆锥角膜、扁平角膜等)做出判别，以及人工晶体植入术之前对植入度数的测定，都需要角膜曲率计测定这些屈光参数。因此，角膜屈光参数的准确性具有特别重要的临床意义。

传统的光机式角膜曲率仪1856年首先由Helmhotlz设计［1］，它是一种基于目视手动调焦原理的用传统刻度显示结果的仪器，主要使用短距望远镜来观测角膜反射的像，测量曲率半径［2］，但有测量速度慢、依赖操作经验等缺点，尤其对年龄大、不配合的人群，测量角膜参数特别不方便。近年来，随着科学技术的发展，自动光机电测量式角膜曲率计(下面简称“自动曲率计”)，以其快速、准确和客观，成为国内外研究的热点。20世纪90年代末，日本Topcon公司和Nidek公司相继研制成功了自动曲率计;21世纪初，韩国佳乐普公司也在市场销售这种仪器。但是，都无文献资料记载技术处理方法，而我国只能依靠进口，价格昂贵。自动曲率计由CCD光学成像系统、电子图像处理系统和机械系统组成，笔者根据国家计量检定规程 JJG1011—2006《角膜曲率计》，提出角膜检测原理、电路系统工作原理和图像处理及计算方法。

1 系统检测原理

对人眼解剖研究的结果表明，角膜是非球面的，但中央表面近似是一个凸面镜［2－3］。本系统的检测原理就是利用人眼的这一光学特点，将排列在一固定大小圆环上的红外LED照射其相同直径的圆靶环(形成靶环系统)，产生圆形光环照射眼睛角膜，角膜反射后经光学系统，成像在 CCD上。由于眼角膜的屈光状态不一致，测量圆环在CCD上成像的大小、形状也不一样。如果是单纯近视或远视眼，在CCD上呈大小一定清晰的圆环;而若是散光眼，则呈清晰的椭圆。通过电子系统对其处理，便可快速、准确、客观地自动测量人眼角膜屈光参数。

测量系统如图1所示，圆靶环上面的光源发出的光线照射到人眼角膜表面。由于角膜表面相当于一凸面镜，会在眼球内部呈正立缩小的虚像，该虚像作为成像系统的物，其形状会因人眼角膜曲率半径的不同而不同，该物通过光学系统在CCD面上成倒立的实像。根据CCD上实像的形状相对于原标准靶环的形状对比即可计算人眼角膜曲率半径R。

图1 角膜屈光参数测量光学原理Fig.1 Optical schematic diagram of cornea diopter parameter

根据CCD上实像的尺寸y0和成像系统的放大倍率β，即可得到靶环虚像的尺寸y′，根据虚像的尺寸即可计算人眼的角膜曲率半径，如图2所示。

图2 角膜曲率半径计算示意Fig.2 Schematic diagram calculating radius of cornea

根据几何关系可知

而

因此得到

式中，R为角膜曲率半径，d为靶环到眼角膜顶点的距离，y为靶环的实际尺寸，y0为CCD上像的尺寸，β为成像系统的放大倍率。从上述计算方法可知，角膜曲率与d、y、y0、β有关。实际测量时 y的误差很小，而β为常数，因此主要考虑的因素为 d，以及y0计算的准确性。d可以通过CCD成像系统的清晰度保证其精度。角膜屈率半径的精度是±0.02 mm，对应的CCD光学成像分辨率为0.1像素，y0的计算要求在亚像素精度才能满足计算的要求。y0的准确度计算也是本研究的重点。

式(3)是假设眼睛角膜仅有球镜的情况，这时y0就是成像在CCD上圆环的直径。如果角膜有散光，就要测量互相垂直方向的曲率半径，这时就相当于计算成像在CCD上椭圆环的长轴和短轴。

求出曲率半径，可算出角膜屈光度值为

式中，D为角膜屈光度值，n为角膜折射率(一般为1.337 5)，r为角膜曲率半径。

将椭圆长、短轴代入式(4)就可算出角膜水平方向屈光度D1和角膜垂直方向屈光度 D2，角膜的散光用C表示，即

散光的轴角是 θ，即是椭圆的长轴或短轴角度。

2 电路系统工作原理

图3是电路系统原理框图。首先，当眼睛放在仪器的前面时，调节机械系统，在CRT监视器上出现眼睛。电路在CRT中心产生方框，眼睛对准方框中心说明光学系统与角膜中央同轴，前后调节大圆环清晰就可以测量。然后，电路在 DSP TMS320F2812(其主频150 MHZ)处理系统的控制下，CCD上成的图像经SAA7111 A/D转换为数字信号，通过CPLD、SRAM图像采集系统，完成角膜CCD图像采集、图像数据处理，系统也采用 CPLD、SRAM、同步分离电路，实现数字测量结果的信号输出，与模拟视频角膜监视信号叠加，在医用高清晰CRT上显示结果，监视眼睛测量。另外，系统由恒流源驱动电路控制红外LED，照射靶环保证测量环的亮度一致。在分别测量两个眼睛时，机械系统带动电位器，记录了两个眼睛的位置，通过DSP上的A/D转换成数字信号，计算出人眼瞳距。角膜屈光参数也可以用微型热敏打印机打印出来，或者通过USB接口把数据保存到计算机上。

图3 电路系统原理Fig.3 Schematic drawing of circuit system

3 屈光参数计算方法

根据光学成像原理，系统计算精度高，必须采用亚像素计算方法［4］。亚像素计算方法有插值、最小二乘曲线拟合等方法，采用大面积插值增加了数据量和数据空间，全部数据参与计算增加了计算量，计算繁琐且计算速度慢。对于平面上任意位置无固定中心的椭圆，计算其长轴、短轴和轴角，目前在欧美国家被公认的较高精度的计算方法是最小二乘法。对于最小二乘法，可以用重心法求出中心，再用中心确定最小二乘法计算。这在理论上可行，但由于有眼睫毛遮挡等因素，这种方法求出中心有误差，造成计算结果不准确;或者环内全部点都参与计算，这种方法精度能满足要求，但数据量大，计算速度慢;结合实际人眼反射的角膜膜环质量好的特点，通过实验，笔者提出边缘检测内外环，再用最小二乘法计算，可以满足测量要求。

图像阈值分割的方法很多，笔者做了多种实验，使用局部最大类间方差阈值等图像分割处理技术［5］，降低计算数据量，求取局部平均阈值，再作为全局阈值的方法，对图像进行二值化分割。在二值化分割后，如果直接计算椭圆环长轴和短轴，重复进行几千个数据量计算，在DSP处理上速度是很慢的［6－7］。因此，程序进行二维梯度计算扫描，可检测出椭圆环的轮廓，同时判别两边缘点的距离，把内外椭圆环分成两个区域存储。实际上，为了提高速度，处理圆环仅取区域一定范围的数据即可，二值化分割和二维梯度计算扫描是同时完成的，也不必实际进行梯度的加减运算，通过程序判别即可，如图4和图5所示。然后，用最小二乘曲线拟合计算椭圆长、短轴和轴角。

图4 椭圆环Fig.4 The ellipse ring

图5 椭圆环边缘轮廓Fig.5 The edge of ellipse ring

选用边缘检测的椭圆最小二乘曲线，拟合计算椭圆长、短轴和轴角。椭圆有内环和外环，分别计算内环和外环椭圆长、短轴和轴角，然后求平均值，就得到椭圆最终的参数。

椭圆曲线的一般方程［8］

其最小二乘拟合方程组为

解线性方程组就可求出 A、B、C、D、E的值。由此可求出椭圆长l1、短轴 l2和轴角θ。

然后，求内环和外环的平均值，得到椭圆最终的参数。椭圆长、短轴乘以一个系数k就是角膜的曲率半径，轴角就是角膜散光的角度。

4 实验

4.1 实验方法

4.1.1 研制的测试仪器

本实验样机如图6所示。测量CRT监视画面，中间方框是测量对准区域，测量时把眼睛放在方框内前后对清楚，即可测量。可以对角膜曲率半径、角膜屈光度进行测量，两者通过程序编成，只要按键即可转换。通过对计量院标准角膜眼进行了测量，符合计量规程。

图6 角膜测量样机Fig.6 Principle instrument of cornea measurement

由于CCD采用ICX405，信号采用SAA7111A/D转换，输出像素频率是13.5 MHz，而CCD信号的时钟频率是9.468 5 MHz，CCD的像元大小是 9.8 μm×6.3 μm，因此DSP处理时 x和 y坐标的比例不一致，需要调整为

图7 角膜CRT测量图示Fig.7 CRT measurement of eye cornea

即y坐标乘1.091，这样修正后计算就没有误差。

4.1.2 《JJG1011—2006角膜曲率计》主要技术要求

1)测量范围:至少能满足6.5～9.4 mm测量范围的要求。

2)精度要求:曲率半径≤8.00 mm，误差为 ±0.02 mm;曲率半径≥8.00 mm，误差为±0.03 mm。

3)测试对象及方法。

实验用计量院专用角膜眼曲率半径为6.67、7.94、9.32 mm和人眼进行。实验时，用两种方法分析比对:第一种用本方法计算结果。第二种方法是使用符合计量院标准的光机式角膜计。国家计量规定检查最小误差为±0.02 mm。

4.2 结果

表1是用所提出的计算方法和光机式角膜计测量的比较情况，计算曲率半径5 mm角膜模拟眼为基准，每组数据测量3次的平均值。表中 l1、l2是相互垂直的角膜曲率半径，θ是l1的角度，标准眼1～标准眼3是计量院专用角膜眼，从上述实验可以看出，采用基于边缘轮廓检测和最小二乘拟合的图像处理方法计算的模拟眼角膜曲率半径与实际基本相同，同光机式角膜计测量也在计量误差范围内，证明该方法可行。另外，每组实验进行了9次测量，重复性在±0.03 mm以内，根据随机误差理论，每组数据测量3次以上求平均，结果在±0.02 mm以内。测量速度在1 s以内，而用光机式角膜计测量时间需要数分钟。

表1 两种方法结果比较Tab.1 Results of the two calculating method

5 讨论

笔者所研究的人眼全自动角膜曲率计与传统手动角膜曲率计相比，具有测量速度快、范围宽、性能稳定、重复性好等优点。

本研究的主要创新点有:

1)系统采用了高景深和自动实时成像光学设计方法，克服了传统光机式曲率计聚焦差和靠测量人员手动调焦主观判断曲率半径情况的缺点，可快速客观地得到测量结果。

2)系统使用面阵CCD相机成像和高速DSP数据处理技术，保证实时成像及数据处理的要求。

3)采用基于边缘特征提取和最小二乘拟合的方法，保证了高测量精度和数据处理速度。

从实验可以看出，测量精度 ±0.02 mm，达到了国家角膜曲率半径的计量要求，并且测量速度相对于传统的光机式角膜曲率计提高了30倍以上，大大地提高了测量的精度和速度。

本研究虽然显著提高了人眼曲率半径的测量速度，克服了传统的光机式角膜曲率计测量中由于眼睛疲劳带来的测量精度下降的缺陷。但是，成像质量仍需人工主观判断，存在一定的由于主观因素带来的误差，因此全自动成像式曲率半径测量闭环系统是未来的主要研究方向。

6 结论

本研究从自动测量人眼角膜曲率半径出发，设计了一套圆环测量系统，采用改进的最大类间方差自适应阈值算法进行二值化图像分割，然后检测椭圆环的边缘，最后使用最小二乘拟合计算椭圆长轴和短轴，达到亚像素精度。实验表明，该方法新颖，对大量数据进行处理，达到快速、准确、实用的目的;该方法达到了计算要求，能够精确计算出人眼角膜屈光参数。

［1］Pantazis M.Visual Instrumentation［M］. New York:McGraw2hill，1999:106 －110.

［2］Maloncy PK.Determination of comed image-foming properties from comal topography［J］.Am J Opthalmol，1987，105 － 223.

［3］Blaker JW.Toward an adaptive model of the human eye［J］.J Opt Soc Am，1980，70(2):220－223.

［4］赵俊奇.人眼屈光度客观式测量的图像采集处理研究［J］.光学技术，2002，28(4):293 －295.

［5］王勤，赵玉环，余景池.图像处理在MDT测量镜片屈光度中的应用［J］. 光电工程，2009，36(12):114－118.

［6］鞠颖，王博亮，谢杰镇.基于裂隙灯显微图像的眼前节特征提取的新方法［J］.中国生物医学工程学报，2004，23(3):193－203.

［7］肖松山，范世福，张思祥.基于人眼视觉特性的图像处理技术［J］.中国生物医学工程学报，2001，20(5):441 －444.

［8］刘书桂，李蓬，那永林.基于最小二乘原理的平面任意位置椭圆的评价［J］.计量学报，2002，23(4):245 －247.