本发明提出一种基于睑板腺图像的全自动腺体分割及量化分析方法，包含以下步骤：导入基于红外成像的睑板腺图像；对睑板腺图像进行R通道提取，干扰信息去除；对图像进行对比度增强处理，二值化形态学处理得到眼睑大概区域；对眼睑区域的上下边界进行线拟合，拟合线映射到增强处理后的睑板腺图像，自动提取眼睑区域；对眼睑区域进行二值化，根据倾斜度及粒子面积过滤信息；计算黏连阈值，根据黏连阈值自动提取黏连腺体；对黏连腺体进行腺体自动分割处理；自动标记每条腺体进行腺体量化分析。本发明提出一种基于睑板腺图像的全自动腺体分割及量化分析方法首先是将睑板腺图像通过图像增强及形态学处理等操作自动得到眼睑的大概区域；然后通过线性拟合得到眼睑上下边界后将边界映射至增强后的睑板腺图像，将睑板腺区域自动分割出来；接着将眼睑图像根据倾斜度过渡信息并计算黏连阈值，根据阈值自动提取黏连腺体并自动分割；最后自动标记每条腺体并进行腺体量化分析，整个过程实现全自动处理。
1.一种基于睑板腺图像的腺体提取方法，其特征在于，包括以下步骤：步骤一：导入基于红外成像的睑板腺图像；
步骤四：对去除干扰信息后的睑板腺图像进行对比度增强处理；
步骤五：对步骤四得到的对比度增强的睑板腺图像进行对比度二值化形态学处理，得到眼睑的大概区域图像；
步骤六：对步骤五中得到的眼睑的大概区域图像的上下边界进行线拟合，拟合线映射到步骤四得到的增强处理后的睑板腺图像，提取睑板腺图像中的眼睑区域图像；
步骤七：对步骤六中得到的眼睑区域图像进行二值化处理并过滤掉无效信息，获得仅有腺体的腺体图像。
2.根据权利要求1所述的一种基于睑板腺图像的腺体提取方法，其特征在于，还包括：步骤八：提取步骤七中得到的腺体图像中的黏连腺体；
3.根据权利要求1或2所述的一种基于睑板腺图像的腺体提取方法，其特征在于，所述步骤三具体为，将R通道提取后的睑板腺图像通过Prewitt算子提取眼睑周边的干扰信息，所述干扰信息至少包括睫毛轮廓和高亮反光点；然后用步骤一中导入的睑板腺图像减去提取的干扰信息。
4.根据权利要求1或2所述的一种基于睑板腺图像的腺体提取方法，其特征在于，所述步骤五具体包括如下步骤：
步骤51：将步骤四得到的对比度增强的睑板腺图像进行二值化并反转处理，通过Reject‑Border去除眼睑周边信息；
步骤52：对步骤51得到的睑板腺图像进行凸包处理和腐蚀处理以分离眼睑区域，保留图中最大面积区域后通过膨胀处理还原原本区域大小，得到眼睑的大概区域图像。
5.根据权利要求1或2所述的一种基于睑板腺图像的腺体提取方法，其特征在于，所述步骤六具体包括如下步骤：
步骤61：将眼睑的大概区域图像通过方向分别为从上往下和从下往上、步进为3个像素的线性扫描找到眼睑的大概区域图像的上下边界轮廓点；
步骤62：由于眼睑边界属于不规则且拐点较多的曲线，因此把上下边界轮廓点分别通过B‑spline的方法线性拟合，分别得到平滑的上下边界拟合线；
步骤64：将步骤63中得到的重构后的眼睑区域图像映射到步骤四中的对比度增强的睑板腺图像，获取睑板腺图像中的眼睑区域图像。
6.根据权利要求1或2所述的一种基于睑板腺图像的腺体提取方法，其特征在于，所述步骤七具体包括如下步骤：
步骤71：将步骤六中提取的眼睑区域图像进行二值化，计算二值化后眼睑区域图像中每个单独粒子的倾斜度和面积以及总粒子的平均面积 及标准差SD，其中，倾斜度SA的计算公式：SA＝arctan k，k为粒子骨架化后拟合直线的斜率；
步骤72：过滤掉眼睑区域图像中粒子的倾斜度不在预定范围内的粒子以及面积小于的粒子。
7.根据权利要求2所述的一种基于睑板腺图像的腺体提取方法，其特征在于，所述步骤八具体为：获取每条腺体的水平段数Nh和垂直段数Nv，计算出水平段数的平均值 垂直段数的平均值 水平段数的标准差SDh、竖直段数的标准差SDv，设水平黏连阈值为垂直黏连阈值为 当腺体存在 或 时，判定为该腺体为黏连腺体。
8.根据权利要求2所述的一种基于睑板腺图像的腺体提取方法，其特征在于，所述步骤九具体包括如下步骤：
步骤91：从步骤八中的眼睑区域图像中提取黏连腺体图像，并将提取的黏连腺体图像在眼睑区域图像中去除，并对黏连腺体通过腺体分割法实现腺体的分割；
步骤92：把分割后的腺体重绘到去除黏连腺体后的眼睑区域图像中并以不同的颜色标记每一条腺体。
9.根据权利要求8所述的一种基于睑板腺图像的腺体提取方法，其特征在于，所述腺体分割法具体包括如下步骤：
步骤911:对黏连腺体图像进行腐蚀操作，判定腺体的个数，如腺体个数为1，再继续腐蚀操作；
步骤912：当腺体个数大于1，进行反转并骨架化，得到腺体的分割骨架图像；
步骤913：将黏连腺体图像减去分割骨架图像，实现黏连腺体的分割；
10.一种基于权利要求1‑9任一项所述的腺体提取方法提取的睑板腺图像的腺体量化分析方法，具体包括如下步骤：
步骤a：对需要分析的腺体进行扫描，扫面方向垂直于腺体拟合直线，扫描为步进为3个像素的线性扫描，通过扫描获取腺体在沿着腺体拟合直线方向上不同位置处的腺体宽度di和腺体中心线长度lMN，通过不同位置处的腺体宽度di获取腺体宽度的平均值D；
步骤b：根据不同位置处的腺体宽度di、腺体宽度的平均值D和腺体中心线长度lMN计算腺体的生物学参数，所述生物学参数包括形变系数、迂曲系数、总腺体的显影值，所述形变系数计算公式为：
其中，l′MN为腺体拟合直线长度，Δαi为腺体中心线相隔3个像素两点的切线夹角，Δsi为两点间的弧长；
其中，GREYi为总腺体的平均灰度值，GREY0为眼睑区域内非腺体部分的平均灰度值。
一种基于睑板腺图像的腺体提取方法及量化分析方法
技术领域
[0001] 本发明涉及睑板腺图像处理领域，尤其涉及一种基于睑板腺图像的腺体提取方法及量化分析方法。
背景技术
[0002] 睑板腺埋藏于眼睑之内侧，开口位于睑缘，具有分泌各种脂质成分，形成脂质层，防止过度泪液蒸发的功能，是维持眼表健康和完整的必要条件。睑板腺功能障碍
(Meibomian Gland Dysfunction，MGD)是一种高发病，MGD通常以终末导管阻塞，腺分泌异
常为特征，常导致眼表上皮损伤、慢性睑炎、干眼病等，严重降低生活质量。虽然国内外研究
者对MGD的患病率进行了一系列的研究，但报道的结果却是相差很大，从3.6％到69.3％不
等，这主要是由于缺乏有效和统一的诊断标准。因此，评价睑板腺的变化对临床诊断和治疗
MGD具有重要意义。
[0003] MGD受到了临床的普遍关注，同时在临床检查中，基于Meibography技术的各种眼表综合分析仪也应运而生。Meibography作为一种红外照明的非接触式成像技术，从整个人
眼睑获取睑板腺图像，能更好地帮助眼科医生评估睑板腺体，对腺体萎缩和丢失情况进行
评估，制定标准化的诊断和治疗策略。通过睑板腺图像，可以直接观察腺体的扩张、扭曲、缩
短、萎缩等形态学变化，直观、定性地进行评价。
[0004] 然而，基于睑板腺图像识别眼睑区域的图像处理方法已有很多，但针对腺体的分割方法以及腺体的量化分析方法却很少。如中国专利CN111008988B和CN111652847A所公开
的睑板腺图像处理方法仅适用于识别眼睑轮廓，其轮廓边界并未进一步拟合，眼睑边界凹
凸不平，所得眼睑区域精度并不高。如中国专利CN109087310A，所公开的睑板腺腺体分割方
法基于深度学习自动分割腺体，分割后对腺体量化的参数仅有面积占比。中国专利
CN111145155A所公开的睑板腺腺体分割方法基于灰度值划分识别腺体，根据扭曲值过滤腺
体，但未对黏连的腺体进行进一步分割，且对腺体量化的参数仅有5种常规生物学信息。这
些参数已经是领域内常用的参数，很多商业机已有该参数的量化，对目前MGD的分级评估没
有先进的突破。
[0005] 本发明提出一种全新的，全自动的眼睑识别与腺体分割图像处理方法，对黏连的腺体进一步分割并对腺体量化的生物学参数提出三个具有临床意义与价值的参数，分别是
形变系数、迂曲系数和显影值。从而进一步减轻医师的工作量，更加客观准确的评价睑板腺
的确实和分级，更加直观地表达腺体的扭曲、膨大、萎缩，MGD诊断的重要特异性体征，为市
场上庞大的睑板腺数据库提供准确的腺体量化分析工具。
发明内容
[0006] 本发明的目的在于提供一种基于睑板腺图像的腺体提取方法及量化分析方法，能够自动处理睑板腺图像，特别是眼睑区域的识别以及黏连腺体的自动分割，并且能够更加
准确地计算出腺体的形变系数、迂曲系数和显影值等生物学参数。
[0007] 为实现上述目的，本发明提供了如下技术方案：
[0008] 一种基于睑板腺图像的腺体提取方法，其特征在于，包括以下步骤：
[0009] 步骤一：导入基于红外成像的睑板腺图像；
[0010] 步骤二：对睑板腺图像进行R通道提取；
[0011] 步骤三：对R通道提取候的睑板腺图像去除干扰信息；
[0012] 步骤四：对去除干扰信息后的睑板腺图像进行对比度增强处理；
[0013] 步骤五：对步骤四得到的对比度增强的睑板腺图像进行对比度二值化形态学处理，得到眼睑的大概区域图像；
[0014] 步骤六：对步骤五中得到的眼睑的大概区域图像的上下边界进行线拟合，拟合线映射到步骤四得到的增强处理后的睑板腺图像，提取睑板腺图像中的眼睑区域图像；
[0015] 步骤七：对步骤六中得到的眼睑区域图像进行二值化处理并过滤掉无效信息，获得仅有腺体的腺体图像。
[0016] 优选地，还包括：
[0017] 步骤八：提取步骤七中得到的腺体图像中的黏连腺体；
[0018] 步骤九：对黏连腺体进行分割。
[0019] 优选地，在步骤三具体为，将R通道提取后的睑板腺图像通过Prewitt算子提取眼睑周边的干扰信息，所述干扰信息至少包括睫毛轮廓和高亮反光点；然后用步骤一中导入
的睑板腺图像减去提取的干扰信息。
[0020] 优选地，所述步骤五具体包括如下步骤：
[0021] 步骤51：将步骤四得到的对比度增强的睑板腺图像进行二值化并反转处理，并通过Reject‑Border去除眼睑周边信息；
[0022] 步骤52：通过convex hull算子进行凸包处理和腐蚀处理以分离眼睑区域，保留图中最大面积区域后通过膨胀处理还原原本区域大小，得到眼睑的大概区域图像。
[0023] 优选地，所述步骤六具体包括如下步骤：
[0024] 步骤61：将眼睑的大概区域图像通过方向分别为从上往下和从下往上、步进为3个像素的线性扫描找到眼睑的大概区域图像的上下边界轮廓点；
[0025] 步骤62：由于眼睑边界属于不规则且拐点较多的曲线，因此把上下边界轮廓点分别通过B‑spline的方法线性拟合，分别得到平滑的上下边界拟合线；
[0026] 步骤63：通过上下边界拟合线重构眼睑区域图像；
[0027] 步骤64：将步骤63中得到的重构后的眼睑区域图像映射到步骤四中的对比度增强的睑板腺图像，获取睑板腺图像中的眼睑区域图像。
[0028] 优选地，所述步骤七具体包括如下步骤：
[0029] 步骤71：将步骤六中提取的眼睑区域图像进行二值化，计算二值化后眼睑区域图像中每个单独粒子的倾斜度和面积以及总粒子的平均面积 及标准差SD，其中，倾斜度SA
的计算公式：SA＝arctan k，k为粒子骨架化后拟合直线的斜率；
[0030] 步骤72：过滤掉眼睑区域图像中粒子的倾斜度不在预定范围内的粒子以及面积小于 的粒子。
[0031] 优选地，所述步骤八具体为：获取每条腺体的水平段数Nh和垂直段数Nv，分别计算出水平段数和垂直段数的平均值和标准差，设水平黏连阈值为 垂直黏连阈值为
当腺体存在 或 时，判定为该腺体为黏连腺体。
[0032] 优选地，所述步骤九具体包括如下步骤：
[0033] 步骤91：从步骤八中的眼睑区域图像中提取黏连腺体图像，并将提取的黏连腺体图像在眼睑区域图像中去除，并对黏连腺体通过腺体分割法实现腺体的分割；
[0034] 步骤92：把分割后的腺体重绘到去除黏连腺体后的眼睑区域图像中并以不同的颜色标记每一条腺体。
[0035] 优选地，所述腺体分割法具体包括如下步骤：
[0036] 步骤911:对黏连腺体图像进行腐蚀操作，判定腺体的个数，如腺体个数为1，再继续腐蚀操作；
[0037] 步骤912：当腺体个数大于1，进行反转并骨架化，得到腺体的分割骨架图像；
[0038] 步骤913：将黏连腺体图像减去分割骨架图像，实现黏连腺体的分割；
[0039] 步骤914：对分割后的腺体进行腺体轮廓修复。
[0040] 本发明还提供了一种基于睑板腺图像的腺体量化分析方法，具体包括如下步骤：
[0041] 步骤a：对需要分析的腺体进行扫描，扫面方向垂直于腺体拟合直线，扫描为步进为3个像素的线性扫描，通过扫描获取腺体在沿着腺体拟合直线方向上不同位置处的腺体
宽度di和腺体中心线长度lMN，通过不同位置处的腺体宽度di获取腺体宽度的平均值D；
[0042] 步骤b：根据不同位置处的腺体宽度di、腺体宽度的平均值D和腺体中心线长度lMN计算腺体的生物学参数，所述生物学参数包括形变系数、迂曲系数、总腺体的显影值，
[0043] 所述形变系数计算公式为：
[0044]
[0045] 所述迂曲系数计算公式为：
[0046]
[0047] 其中，Δαi为腺体中心线相隔3个像素两点的切线夹角，Δsi为两点间的弧长；
[0048] 所述显影值计算公式为：
[0049]
[0050] 其中，GREYi为总腺体的平均灰度值，GREY0为眼睑区域内非腺体部分的平均灰度值。
[0051] 与现有技术相比，本发明具有以下有益效果：
[0052] 1)本发明能够实现对睑板腺图像中的腺体进行自动提取；并且能够对黏连腺体进行分割；2)本发明能够对腺体的各个参数进行分析。
附图说明
[0053] 图1为全自动腺体分割及量化分析方法的示意图；
[0054] 图2为眼睑区域自动识别的图像处理方法示意图；
[0055] 图3为腺体自动提取的示意图；
[0056] 图4为腺体自动分割算法的示意图；
[0057] 图5为单个睑板腺生物参数计算要素示意图
[0058] 图6为代表性的睑板图像全自动腺体分割及量化分析结果示意图
具体实施方式
[0059] 以下结合附图和具体实施方法对本发明做进一步的描述：
[0060] 实施例一
[0061] 如图1所示，本发明提出了一种基于睑板腺图像的腺体提取方法：包括以下步骤：
[0062] 步骤一：导入基于红外成像的睑板腺图像。
[0063] 步骤二：对睑板腺图像进行R通道提取。
[0064] 步骤三：对R通道提取候的睑板腺图像去除干扰信息。具体地，将R通道提取后的睑板腺图像通过Prewitt算子提取眼睑周边的睫毛轮廓和高亮反光点等干扰信息，然后用步
骤一中导入的睑板腺图像减去提取的干扰信息实现去除，如图2b所示。
[0065] 步骤四：对去除干扰后的睑板腺图像进行对比度增强处理。具体地，将干扰信息去除的睑板腺图像通过中值滤波器、高亮细节滤波器和拉普拉斯滤波器组合处理后，得到对
比度增强的睑板腺图像，如图2c所示。
[0066] 步骤五：对步骤四得到的对比度增强的睑板腺图像进行对比度二值化形态学处理得到眼睑大概区域。具体采用如下步骤：
[0067] 步骤51：将步骤四得到的对比度增强的睑板腺图像进行二值化并反转处理，并通过Reject‑Border去除眼睑周边信息，如图2d所示；
[0068] 步骤52：通过convex hull算子进行凸包处理和腐蚀处理以分离眼睑区域，保留分离得到的眼睑区域图像中的面积最大的图像，并将面积最大的图像通过膨胀处理还原成分
离之前的眼睑区域的大小，得到眼睑的大概区域图像，如图2e所示。通过分离能够将眼睑区
域边缘的一些不整齐的部分分离，这些分离图像中，中间部分的面积最大，边缘的部分的面
积很小，中间部分最能代表眼睑的图像，去除掉边缘部分，再将中间部分进行膨胀，即可得
到眼睑的大概区域图像。
[0069] 步骤六：对步骤五中得到的眼睑的大概区域图像的上下边界进行线拟合，拟合线映射到增强处理后的睑板腺图像，以提取眼睑区域图像。具体采用如下步骤：
[0070] 步骤61：将眼睑的大概区域图像通过方向分别为从上往下和从下往上、步进为3个像素的线性扫描找到眼睑的大概区域图像的上下边界轮廓点；
[0071] 步骤62：由于眼睑边界属于不规则且拐点较多的曲线，因此把上下边界轮廓点分别通过B‑spline的方法线性拟合，分别得到平滑的上下边界拟合线；
[0072] 步骤63：通过上下边界拟合线重构眼睑区域图像(如图2f所示)，即得到以两条拟合线为边界的眼睑区域图像；
[0073] 步骤64：将步骤63中得到的重构后的眼睑区域图像映射到步骤四中的对比度增强的睑板腺图像，获取睑板腺图像中的眼睑区域图像，如图2g所示。当然，为了更好展示拟合
边界的准确性，也可以把拟合边界映射到步骤一中的睑板腺图像，如图2h所示。
[0074] 步骤七：对眼睑区域图像进行二值化，并根据步骤六中眼睑区域图像中的粒子倾斜度及粒子面积对二值化后的眼睑区域图像进行过滤，所述粒子可以理解为，图像是由像
素组成的，在一个图像中，有些位置像素会聚集在一起形成像素团，这些像素团就是一个一
个的粒子。所述步骤七具体包括如下步骤：
[0075] 步骤71：将步骤六中提取的眼睑区域图像进行二值化，计算图中每个单独粒子的倾斜度和面积以及总粒子的平均面积 及标准差(SD)，其中，倾斜度(Slant Angle)的计
算公式：SA＝arctank，k为粒子骨架化后拟合直线的斜率，所述骨架化为现有算法；
[0076] 步骤72：过滤倾斜度不在区间[45,135]范围内的粒子(如图3b所示)以及面积小于的粒子，得到过滤后的眼睑区域图像，如图3c所示。
[0077] 步骤八：提取黏连腺体，所述黏连腺体为相邻的腺体黏连在一起，需要将黏连腺体分开。具体地，获取每条腺体的水平段数(number of horizontal segments，Nh)和垂直段
数(number of vertical segments，Nv)，分别计算出水平段数和垂直段数的平均值和标准
差，设黏连阈值为 当腺体存在 或 时，判定为该腺体
为黏连腺体，并将黏连腺体标记为与其它腺体不同的颜色，具体如图3c所示。
[0078] 步骤九：对黏连腺体进行分割。具体包括如下步骤：
[0079] 步骤91：提取黏连腺体图像(如图4a所示)，并将提取的黏连腺体图像在眼睑区域图像中去除(如图3d所示)，黏连腺体通过腺体分割法实现腺体的分割；
[0080] 步骤92：把分割后的腺体重绘到去除黏连腺体后的眼睑区域图像中并以不同的颜色标记每一条腺体，即得到了分割后的腺体，如图3e所示，图3e中实际为彩色。
[0081] 所述腺体分割法具体包括如下步骤：
[0082] 步骤911:对黏连腺体图像进行腐蚀操作，判定腺体的个数，如腺体个数为1，再继续腐蚀操作；
[0083] 步骤912：当腺体个数大于1(如图4b所示)，进行反转并骨架化，得到腺体的分割骨架图像(如图4c所示)，所述反转为将图像中的0变为1,1变为0；
[0084] 步骤913：将黏连腺体图像减去分割骨架图像，实现黏连腺体的分割；
[0085] 步骤914：对分割后的腺体进行腺体轮廓修复以修复腺体部分区域受分割骨架影响的形成的缺损(如图4e所示)。具体地，如图4f所示，分别对分割后的腺体进行外轮廓提取
并填充，最后减去腺体的内轮廓得到修复后的腺体，所述内轮廓的获取采用现有算法。
[0086] 如图4h所示为腺体自动分割算法的效果图。
[0087] 实施例二
[0088] 该实施例为对腺体的量化分析方法。
[0089] 步骤a：对需要分析的腺体进行扫描，扫面方向垂直于腺体拟合直线，扫描为步进为3个像素的线性扫描，通过扫描获取腺体在沿着腺体拟合直线方向上不同位置处的腺体
宽度di，通过不同位置处的腺体宽度di获取腺体宽度的平均值D；
[0090] 步骤b：根据不同位置处的腺体宽度di和腺体宽度的平均值D计算腺体的生物学参数，包括长度，直径，面积，形变系数和迂曲系数，并计算总腺体占眼睑的面积比，生物学参
数均值和总腺体的显影值。如图5所示为单个睑板腺生物参数计算要素。
[0091] 进一步地，所述形变系数(Deformation Index，DI)计算公式为：
[0092]
[0093] 其中，di为腺体线性扫描方向上两边的距离，D为腺体两边平均距离；
[0094] 进一步地，所述迂曲系数(Tortuosity Index，TI)计算公式为：
[0095]
[0096] 其中，lMN为腺体中心线长度，l′MN为腺体拟合直线长度，Δαi为腺体中心线相隔3个像素两点的切线夹角，Δsi为两点间的弧长；
[0097] 进一步地，所述显影值(Signal Index，SI)计算公式为：
[0098]
[0099] 其中，GREYi为总腺体的平均灰度值，GREY0为眼睑区域内非腺体部分的平均灰度值。
[0100] 最后，睑板腺图像以及分析的结果按模板捆绑保存，通过数据导出以Excel，png，bmp，jpeg，txt和PDF等多个常用的数据格式导出数据，为临床医生提供可靠的数据，方便临
床医生的随访研究。图6展示的是代表性的睑板图像全自动腺体分割及量化分析结果。
[0101] 显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可
以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本
发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求
的保护范围之内。
