一种角膜Placido圆环边缘检测方法

郭雁文

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

0 引言

角膜地形图参数的精确计算，对于眼科医学应用具有重要的意义，而角膜地形图预处理的精度将决定参数的计算精度［1-2］.在角膜地形图的获取过程中，由于受光斑、睫毛等不同噪声的特异性的影响，就会导致在直角坐标系中无法获得较高精度的预处理结果.对临床诊断来说，角膜地形图的处理有三个重要的因素［3］:①角膜地形图重建的精度;②角膜覆盖度(表示为最大角膜直径)理想状况应该是覆盖整个角膜(即眼球中黑色虹膜的部分);③对于眼角膜实际拍摄捕捉过程中存在的各种问题的解决能力，例如被检测者头部的晃动，眼球的移动，泪膜的质量，粘液的存在，特别是睫毛的遮挡等造成的影响［4-5］.而要有效地解决这些问题，首先应先对CCD 相机捕捉到的角膜圆环图像进行预处理，从而找到图像重建方案［6］.

现有的基于直角坐标系的预处理算法除计算量较大之外，还会导致两个极端:①导致有用的数据缺失严重;②人工拟合的数据过多，在进行参数计算前无法保证数据的精度［7-8］.因此，根据角膜地形图原始图像的特性，本文提出了一种基于极坐标的原始图像的处理算法，在图像坐标转换完成后，重点研究在极坐标系下角膜地形图的自适应边缘检测及边缘生长，有效地克服了原始图像的数据丢失和获取的数据不完全而导致的误差问题.

1 传统图像预处理算法存在的问题

图1 为实验中所获取的基于Placido 环的角膜地形原始图像.由于患者眼睑很难张开到所要求的宽度，泪膜不稳定，角膜表面的可涵盖粘液等，都导致了所获取的图像信息丢失严重［9-10］.尤其是睫毛，大量遮盖了有用的信息［11-12］.如图1 所示，眼睑难以张开，睫毛又浓厚，这些都会导致所捕捉到的角膜地形图质量很差.如果利用传统的基于直角坐标系的预处理算法，那么对其直接进行预处理所得到的圆环图就会受到睫毛等噪声信号的严重影响，难以进行下一步的参数计算，如图2 所示.

根据直角坐标和极坐标系的特点，直角坐标系中的一个圆，如果以圆心作为极点，转换到极坐标系后将是一条直线［13］.从图1中可以看到，Placido 环图像呈现为准圆形的同心圆，因此如果能找到这些同心圆的圆心，将其作为极点，将Placido 环图像转换为极坐标域，那么这些亮环将成为单一方向的准直线型，这样对于去除眼睑与眼皮的工作，以及有效检测亮环的边缘工作，都将会被大大地简化.基于此，本文提出了一种新的基于极坐标的角膜图像预处理算法.

图1 Placido 原始图像Fig.1 Original image of Placido

图2 传统算法预处理结果Fig.2 Traditional algorithm preprocessing results

2 基于极坐标系的角膜地形图像自适应预处理算法

基于极坐标系的角膜图像预处理，首先要完成坐标系的转换，同时对转换坐标系之后的图像进行重新生成，然后在新坐标系环境下对图像进行进一步处理.改进算法的具体流程如图3 所示.

图3 新算法工作步骤示意图Fig.3 Steps for implementation of the new algorithm

2.1 极坐标图像生成

在实际拍摄中的Placido 图像，图像的中心往往不是圆环的中心，因此圆环中心的定位非常重要.利用原始图像的准圆环特性，设原角膜地形图像为f(x，y)，nx，ny分别是该图像在X 方向与Y 方向的像素点数.定义一个方形子图像fES(m，n)，对子图像fES(m，n)实现Hough 变换检测圆心，利用原始图像近似同心圆的特征，当检测到两个不同半径的圆，并且其圆心相同时，就可以确定圆心，子图像圆心的定位结果如图4 所示.利用该圆心作为新极坐标系的极坐标原点，将原始图像重新调整和生成为以极坐标原点为中心的方形新图像，目的是为了在下一步进行极坐标转换后，圆环能被有效地展开为直线的形式.矫正后的图像如图5 所示.

图4 子图像圆心定位结果Fig.4 The results of sub image central location

图5 定位后重新生成的图像转换至极坐标系后如图6 所示.图6 中下方如山丘一样的黑色弧形为重定位时所扩展的部分.

图5 矫正后的图像Fig.5 The corrected image

图6 转换极坐标后的图像Fig.6 Image of polar coordinate conversion

2.2 基于极坐标系的自适应水平边缘检测

从图6 中可以看到，有用的圆环呈现近似横线的形式，而睫毛等噪声信号以纵向线条或不同的角度存在，因此复杂的噪声干扰去除就变为提取水平方向的梯度算子，从而大大降低了噪声的影响.因此，先对极坐标转换后的角膜地形图像进行水平边缘检测.自适应水平边缘检测步骤如图7 所示.

图7 自适应水平边缘检测步骤Fig.7 Steps of self-adaptive horizontal edge detection

首先对图像进行平滑处理.选择自适应维纳滤波器，使原始图像和其滤波后的图像之间的均方误差达到最小，根据局部方差来调整滤波器效果.利用Sobel 算子对图像进行水平方向的检测.非极大值抑制实质上是找到边缘检测强度数据中的最高点，利用上述水平边缘检测的结果，对边缘方向法线两侧的两个点插入边缘强度［14］，如果感兴趣点的边缘强度比这两个值大，则返回该值，否则丢弃.当非极大值抑制定义的边缘点大于上限阈值时(或小于上限阈值时)，标记为边缘点，然后利用8 连通域查找它的相邻点以确定它们是否也大于上限阈值(或小于上限阈值)，如果形成相邻点连起来的支径，则为检测到的边缘图像，处理结果如图8 所示.

图8 自适应水平边缘检测结果Fig.8 Results on self-adaptive horizontal edge detection

2.3 形态学后续处理

如图8 所示，在角膜地形边缘图像中，因为睫毛、光斑等影响而丢失了部分有用的信息，导致图像中线段有短缺区域.而在形态学算子中，利用结构元素对图像做膨胀运算，可以填充目标内部狭窄的裂缝和长细的窄沟，消去小的孔洞，而腐蚀运算可以消除图像中小的成分［15］.这非常适合用来对原始角膜地形图像进行有效的弥补工作与噪声消除，最后对圆环周围残余的睫毛及眼睑进行连通区域处理，转换回直角坐标系中即可.

如果选择结构元素SE1(x，y)为线性算子，那么膨胀可填充图像中的细长窄沟，对极坐标下的角膜地形图中的准直线具有弥补作用.而腐蚀可以消除图像中小的成分，并将图像缩小，从而使其补集扩大.因此可选用菱形结构元素SE2(x，y)对图像进行残余噪声去除.对极坐标边缘图像fE作膨胀运算，可表示为

最后标记连通分量，计算连通面积，面积过小的像素点去除掉.这里，为了避免损坏有用的信息，对于连通区域的查找采用除水平方向的连通区域查找，结果如图9 所示.由图9 可以看到，因为睫毛影响所造成的数据损伤已基本修复，大的数据短缺不应全部补齐，否则会影响到数据精度，最后将图像转换到直角坐标系中.

图9 形态学后续处理结果Fig.9 Results on morphological follow-up processing

图10 为处理结果的比较，其中图10(a)为原始图像矫正后的图像，图10(b)为图10(a)红线处的灰度值.可以看到，由于光斑的影响，图10(b)中有多处圆环漏检.图10(c)为处理结果图像，图10(d)为相同行的边缘检测结果图，可以看到圆环均能有效地检测出来，这充分说明了自适应预处理算法使得数据精度有了大幅度提高.

图10 最终处理结果比较Fig.10 Compare of final processing results

3 实验验证

图11 为实验中获取的三组被测者角膜地形图像及处理结果显示，其中原始图像有光斑影响，中心点错位，睫毛影响等共同存在的现象.

经过自适应预处理后，从处理结果可以看到，图(a)中被测者甲的图像受睫毛的影响较为严重，导致直接获取边缘图像时数据损失严重;图(b)的处理结果能够去除其大量内外睫毛的影响，实现了如图(c)中三维显示的较为完整的角膜建模图像.图(d)为被测者乙的角膜原始图像，其中心点严重错位，且有光斑影响，直接进行参数计算误差很大;图(e)能准确矫正圆心位置，并准确获取其有效数据.图(g)中被测者丙的角膜原始图像光斑影响严重，直接检测时圆环图像有严重的断裂现象;图(h)中有效去除了光斑的影响，弥补了由光斑导致的数据确失的现象，检测结果数据损失小，且精度较高，能满足后续测量参数计算的要求.

图11 三组原始实验数据及其处理结果Fig.11 Three sets of original experimental data and their processing results

4 结论

本文根据Placido 角膜原始图像的特点，采用极坐标图像的转换方式，把复杂的噪声干扰去除方法转化为极坐标系下的水平边缘检测和边缘跟踪问题，大大简化了角膜地形图原始图像的预处理难度，同时增加了算法的处理精度，获得了圆环边缘检测完整的结果图像，满足了后续参数计算的要求.大量实验数据证明:该算法具有数据损失小、精度高，算法易实现的特点，满足眼科医学的实际使用要求.

［1］刘祖国.角膜地形图学［M］.广州:广东科技出版社，2001.

［2］Alkhaldi W，Iskander D R，Zoubir M.Model-order selection in Zernike polynomial expansion of corneal surfaces using the efficient detection criterion［J］.IEEE Transactions on Biomedical Engineering，2010，57(10):2429-2438.

［3］Zhu Mingxia，Collins M J，Yeo A C.Stability of corneal topography and wavefront aberrations in young Singaporeans［J］.Clinical and Experimental Optometry，2013，96(5):486-493.

［4］Wizert A，Iskander D R，Cwiklik L.Organization of lipids in the tear film:a molecular-level view［J］.Journal of Biomedical Optics，2014，9(3):e92461.

［5］Atwood D K，Andersen H E，Matthiss B，et al.Impact of topographic correction on estimation of aboveground boreal biomass using multi-temporal，l-band backscatter［J］.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2014，7(8):3262-3273.

［6］Alonso-Caneiro D，Iskander D R，Collins M J.Estimating corneal surface topography in videokeratoscopy in the presence of strong signal interference［J］.IEEE Transactions on Biomedical Engineering，2008，55(10):2381-2387.

［7］Polette A，Auvinet E，Mari J L，et al.Construction of a mean surface for the variability study of the cornea［J］.IEEE Conference Publications，2014:328-335.

［8］周洪亚，沈建新，高绍雷，等.基于精确中心定位的角膜复原算法［J］.生物医学工程学杂志，2011，28(5):872-875.Zhou Hongya，Shen Jianxin，Gao Shaolei，et al.An algorithm of corneal reconstruction based on precise location of corneal center［J］.Journal of Biomedical Engineering，2011，28(5):872-875.(in Chinese)

［9］Snellenburg J J，Braaf B，Hermans E A，et al.A forward ray tracing for image projection prediction and surface reconstruction in the evaluation of corneal topography systems［J］.Optics Express，2010，18 (18):19324-19338.

［10］Florindo J B，Soaresa S H M，de Carvalho L A V，et al.Mumford-Shah algorithm applied to videokeratography image processing and consequences to refractive power values［J］.Computer Methods and Programs in Biomedicine，2007，87(1):61-67.

［11］Carvalho L A.Absolute accuracy of placido-based videokeratographs to measure the optical aberrations of the cornea［J］.Optometry and Vision Science，2004，81(8):616-628.

［12］Alkhaldi W，Iskander D R，Zoubir A M，et al.Enhancing the standard operating range of a placido disk videokeratoscope for corneal surface estimation［J］.IEEE Transactions on Biomedical Engineering，2009，56，(3):800-809.

［13］Thalji Z，Alsmadi M.Iris recognition using robust algorithm for eyelid，eyelash and shadow avoiding［J］.World Applied Sciences Journal，2013，25(6):858-865.

［14］Nixon M S.特征提取与图像处理［M］.第2 版.北京:电子工业出版社，2010:109-115.

［15］Burger W，Burge M J.Principles of Digital Image Processing:Advanced Methods (Undergraduate Topics in Computer Science)［M］.New York:Springer，2013.