骆仲舟，邓宇晴，王耿媛，黄远聪，林伟锋，袁进

(中山大学中山眼科中心工程中心，广州 510000)

睑板腺是人体最大的皮脂腺，埋藏于上下睑板之中，其开口位于睑缘，从眼睑边缘垂直穿过后眼睑至睑板的相对边缘。睑板腺的功能是合成、储存并分泌多种脂质，所分泌的脂质能构成泪膜的脂质层，直接影响泪液蒸发及泪膜的稳定性，具有重要的生理功能[1-2]。睑板腺的功能与其形态相关，炎症、堵塞、激素水平变化、环境因素等多种原因均可以影响睑板腺的数量和形态，进而改变睑板腺油脂分泌功能，导致泪膜稳定性下降，发生睑板腺功能障碍(meibomian gland disease，MGD)。

MGD是蒸发过强型干眼症的最常见的病因[3]。蒸发过强型干眼症约占据干眼构成的86%，严重时会导致严重和慢性的症状以及继发性结膜和角膜的改变，是一种困扰人们视觉质量乃至生活质量的常见病[4]。因此，睑板腺腺体的扭曲、膨大、萎缩乃至缺失是MGD诊断的重要特异性体征[5]。

目前，在临床上通过睑板腺红外成像技术使得眼科医生能直接观察到睑板腺腺体结构，可对腺体萎缩和丢失情况进行粗略评估[6]。但如何对上述腺体形态学变化进行量化分析，值得深入研究。2012年，Srinivasan等[7]发表了一种利用ImageJ Software (National Institute of Health；http://imagej.nih.gov/ij)计算MGD患者的睑板腺缺失率的方法，该方法是用ImageJ软件的自由选择工具勾画出睑板腺缺失的面积进行测量，计算MGD患者睑板腺缺失的面积与整个睑板腺的面积的关系，并用百分比表示睑板腺缺失率。然而，使用该软件用户必须手动定义每个患者腺体区域。不同的操作者可能会画出不同的腺体区域，导致结果存在一定的差异。2013年，Arita等[8]对睑板腺图片自动确定测量区域，并对上眼睑各睑板腺的轮廓进行分析。这种方法由于只分析睑板腺本身的轮廓，而没有计算睑板腺的面积[8]，虽对分析睑板腺细微变化有利，但没能客观的表现睑板腺缺失的程度。在针对基于红外成像技术的睑板腺图片和睑板腺的生物参数，既要准确分析睑板腺腺体的形态生物参数和精确计算睑板腺缺失的百分比，本文章自主研发了睑板腺腺体自动分析软件(Meibomian Gland Bio-image Analyzer V3)，并首次将该技术用于计算睑板腺腺体长度、宽窄、面积百分、形变系数及显影值等生物参数。精确的睑板腺量化分析系统，可有效帮助MGD的临床诊断及分级。

1 材料与方法

1.1 基于红外成像原理的睑板腺图像

本项目使用的红外成像设备是Oculus Keratograph 5M，是德国Oculus公司的眼表综合分析仪设备，其睑板腺拍摄功能是通过红外透射睑板腺拍摄，具有增强对比模式，重点突出腺体，更易于辨认。其设备采用的是数码CCD照相机，其分辨率为1 360×1 024像素，经过增强对比模式后，高清睑板腺图像分辨率为1 088×512像素。工作距离100 mm，使用中心波长为830 nm的红外光源。使用低反射栅格测试板(Edmund #62-949，2×2英寸，0.0625 mm点直径，0.125 mm点间距)测量出，1个像素表示实际尺寸为0.025 mm ×0.025 mm。

1.2 睑板腺图像分析系统

睑板腺图像量化分析分析系统是基于LabVIEW 2016(National Instruments Corporation，US)和NI视觉模块开发，操作系统为Windows 7及更高版本。系统分为图像识别，图像预处理和生物参数分析3大模块。图1所示为一例完整的睑板腺图像量化分析系统界面图。

1.2.1 睑板腺图像识别

由于Oculus Keratograph 5M(以下简称为K5)所导出图像的方式是通过截取屏幕，每个设备使用的分辨率，软件版本并非统一，因此增强对比模式的睑板腺图像在导出的图像上位置不是固定的。针对K5所导出的图像(图2)，系统需要自动识别增强对比模式的睑板腺图像。根据导出图像的灰度差和图片位置，可以发现所需要识别的增强对比模式的睑板腺图像相对边框较暗(灰度值较低)且位置在图像中心点的左上方。因此第1步需要对图像进行分块处理；第2步根据位置的选择需要识别的区域，通过Graham算法[9]截取有效图像(图2D)。

1.2.2 睑板腺图像预处理

使用K 5 的增强对比模式来摄取图像，自动把采集的睑板腺图像增强了睑板腺和眼睑的对比度，通过识别得到分辨率为1 088×512像素的高清睑板腺图像。但K5设备虽把睑板腺和眼睑的对比度增强，更容易观测，但仍然无法实现自动分割，有效提取睑板腺腺体。因此，需要对睑板腺图像进行图像预处理，自动提取出腺体轮廓。图像预处理核心算法部分是基于图像算法Convolution和Morphology算法改写，结合形态滤波及Convex Hull算法，其实现步骤如下。

步骤一：由于睑板腺和眼睑的平均灰度为190和140，而眼睫毛和高亮发光点平均灰度为20和250，因此利用灰度差，通过控制二值化的阈值和边缘检测，快速识别睑板腺的眼睫毛，高亮反光点和睑板腺的轮廓，再使用形态学滤波，可以去除这两部分的图像信息。

其中形态滤波方式如下：A为原图矩阵，B为结构矩阵，R为处理后矩阵。

步骤二：去除眼睑毛和高反光点后，容易精准识别出睑板腺的轮廓，通过对Convolution和Morphology改写算法，结合形态学滤波和Convex Hull算法[10]，去除无效边缘信息，得到仅剩睑板腺的图像。

改写算法核心如下，f为原图矩阵，b为结构矩阵，h为核矩阵，G为处理后矩阵：

图1 睑板腺图像量化分析系统界面图Figure 1 The interface of Meibomian gland image quantitative analysis system

图2 图像分块与识别图像处理图Figure 2 The figure of the image segmentation and recognition image processing

其中公式[2]、公式[3]中，Df和Db分别表示图像输入矩阵f(x,y)和结构矩阵b(x,y)的定义域。公式[4]中，Dh表示核矩阵h(x,y)的定义域。

步骤三：得到睑板腺图像后，进行图像处理重建，增强腺体的对比度，分割腺体与眼睑，自动提取腺体形态。其中f为原图矩阵，k为结构矩阵，M为处理后矩阵；

步骤效果图如图3所示。

通过上述的方法可以在图像预处理的步骤中实现睑板腺区域与腺体的自动提取。但考虑到实验中仍有部分图片无法自动提取，软件同时也开放了手动勾画的模式，操作者可以使用手动勾画的方式把睑板腺的区域勾画出来，系统再根据所勾画的区域提取出睑板腺区域和腺体。

1.2.3 睑板腺的生物参数

睑板腺形态与其功能相关，炎症、堵塞、激素水平变化、环境因素等多种原因均可以影响睑板腺的数量和形态，进而改变睑板腺油脂分泌功能，导致泪膜稳定性下降，发生MGD。因此睑板腺腺体扭曲、膨大、萎缩乃至缺失是MGD诊断的重要特异性体征。

根据腺体的形态，可以计算腺体的长度、直径、面积、形变系数(deformation coefficient，DC)、显影值和中央腺体的腺体百分比6个生物参数。通过检测腺体的边缘，从上到下，每4个像素1个步进计算扫描腺体，计算腺体两边界点A，B的坐标。通过对点A，B的计算，计算出生物参数。计算公式如下：

图3 睑板腺图片预处理过程图Figure 3 The figure of the meibomian gland image preprocessing process

其中，x，y分别表示像素矩阵的横纵坐标，i为二值化像素值。

睑板腺异常轻中度患者，腺体普遍存在扭曲和变形，形变系数是腺体的形体参数。根据临床上腺体的病理性改变，基于弧弦比模型，建立形变系数生物参数。腺体的变异程度量化分析可进一步指导患者病情分级，同时可为进行个体精准化治疗提供评价标准。

其中P为腺体单边周长：

通过对形变系数的比较(图4)，越趋于正常的腺体，形变系数越小，腺体的扭曲和堵塞都会影响形变系数。

从解剖结构上看，睑板腺由腺体导管和腺泡组成，大大小小的腺泡依附在腺体导管的周围，腺泡不停地分泌一种油脂和蛋白的混合物，通过腺体导管不断往外排出，最终到达眼睑的开口处[11]。混合物随着每一次的眨眼动作，在泪液表面自由扩散，涂抹在泪膜的表面，形成脂质层。作为泪膜表面的活性剂，构成一个泪膜屏障，使泪膜的蒸发减慢。当脂质层的稳定性下降，泪膜蒸发加快，维持时间短，可导致蒸发过强型干眼症。因此，通过显影值来表示睑板腺中腺泡的含量。

在同一的照明条件下，不同组织吸收光线的数值叫做光密度(optical density，OD)。利用摄像机采集光信号通过光电转换原理，在图像中以灰度值来表示，测量目标吸收的光线越多，光密度数值越大[12]。当光线照射一个密度均匀的物质时，被吸收光能量的多少与该物质的含量成正比，而与照射光的强度无关[13]。因此，光密度可以用来衡量睑板腺中脂质的含量。以光密度为模型，定义睑板腺的显影值，表示腺体的脂质含量。

显影值公示如下：

图4 腺体的生物参数模型与分析对照图Figure 4 Comparison diagram of the biological parameter model and analysis of the gland

图5 腺体的显影值分析对照图Figure 5 Contrast figure for development of the glands

在临床工作中，评估患者的腺体情况以中央5～8条腺体为主，且睑板腺检查中以中央分腺体最为清晰，中央腺体百分比及长度评估结果可代表患者主要腺体的功能情况，直径可部分反映腺体的阻塞程度。

1.3 统计学处理

采用SPSS 13.0统计软件包进行数据统计分析。计量资料以均数±标准差()表示。重复性评估为同一图像来源不同操作者的2次测量值，采用重复观察平均差异()、变异系数(coefficient of variation，CV%)及组内相关系数(intraclass corr elation coefficient，ICC)表示，CV%＜10%时重复性较好，统计学上认为ICC＞0.75时具有良好的复测可信度。两次分析结果差异用配对样本t检验。P＜0.05为差异有统计学意义。

2 结果

2.1 临床实验设计

由中山眼科中心招募正常受试者24 人(42 眼)，年龄(34±8)岁，其中，男性受试者10人，女性受试者14人。采用Oculus Keratograph 5M进行睑板腺红外摄影。其中选取受试者中的10人(20眼)在同型号的设备上由两名已接受培训的操作员分别行睑板腺红外摄影，所拍摄的图像数据通过睑板腺分析系统(Meibomian Gland Bioimage Analyzer V3)进行生物参数分析。本实验设计手动选择上睑结膜中央5条腺体进行形态学参数的定量分析，对有2名操作员所拍摄的数据进行重复性测试，分析其数据在眼表疾病中的诊断价值。

2.2 实验结果

通过对24位受试者(42眼)的睑板腺图像进行分析，睑板腺生物参数测量均值：腺体平均直径为(0.48±0.09) mm，腺体平均长度为(5.25±0.68) mm，腺体平均面积为(2.12±0.53) mm，腺体形变系数分别为10.01±3.85，显影值为6.32±1.23，中央5条腺体占中央区域面积百分比为(10.94±2.20)%，腺体占上睑结膜面积百分比为(58.07±8.13)%。

选取受试者中的10 人(20眼)在同型号的设备上由第2 名操作员进行第2 次睑板腺红外摄影进行重复性的校验评估，通过配对样本t检验的方法得出结果。所有指标的P值均大于0.05，CV均小于5%，ICC均大于0.95。两次分析的各生物参数差异均无统计学意义，重复性较好，可信度高。

3 讨论

本文介绍了一种高效的、客观的量化分析软件及睑板腺综合分析系统，可用于快速识别睑板腺的区域并自动提取睑板腺的红外成像中的腺体，自动计算中央每条腺体的长度、直径、面积、形变系数和显影值，并提供中央区域腺体的百分比。还可以手动勾选特定腺体区域进行定位随访追踪，为临床提供了多个生物参数，更客观、更高效、更精准地展现睑板腺形态学的变化以及睑板腺缺失的情况。本研究表明：随着腺体的迂曲和堵塞的变化，变异系数有明显差异，腺体越趋于正常，变异系数越小。在同一光源照射下，显影值高为脂质含量高，分泌脂质功能好；而显影值低为脂质含量低，分泌脂质功能差。显影值不受睑板腺图像的整体光亮度影响，但会受到局部反光的影响。

在实验过程中，由于部分测试者对翻眼睑的不适用，睑板腺图像位置不对，反光不均匀和眼睑暴露不完整，这会导致睑板腺图像导入系统后会提取不精准，小部分的干扰信息无法去除。另外当测试者睑板腺严重缺失的，由于缺失的大区域在光源照射下，偶尔会有不同程度的高反光，这会导致睑板腺图像导入后，系统会判断为腺体，属于无法去除的干扰项。以上两点还需要通过对图像算法的优化和提高，进一步提高系统对干扰项的处理能力，增加系统的精确度。

目前中国有近3亿干眼患者，临床检查需求量巨大，将睑板腺综合分析系统植入到干眼筛查检查项目中，有望为临床上对睑板腺腺体功能评估提供新的非侵入性参考指标，进一步帮助医生进行标准化的诊断及疗效评估。