CN104185858A - 睑板腺成像的计算方法和装置 - Google Patents

Abstract

包含睑板腺的区域的眼部图像被自动处理。通过在眼部图像中使用来获得表征眼部图像中的睑板腺的一个或多个数值参数，使用一个或多个数值参数获得分级，这个过程可以获得指示睑板腺健康的分级。这些数值参数包括一个表征由尺度空间变换得到的重要特征图的尺度参数间差异的参数，和/或通过测量眼部图像中代表各自腺体的线获得的参数。使用伽柏滤波器作为一个局部滤波技术，睑板腺可以在眼部图像中被识别。伽柏滤波器的形状、局部空间支持和取向的参数化对于检测睑板腺都很重要。

Description

睑板腺成像的计算方法和装置

技术领域

本发明涉及一种用于处理睑板腺图像来获得表征腺体异常的信息的计算方法和装置，所述图像指示医疗状况。

背景技术

睑板腺为位于睑板内的眼皮边缘处的皮脂腺，用于供应睑脂，睑脂是一种可防止眼睛泪膜的蒸发的油性物质。睑脂是一种防止眼泪溢出到脸颊的脂类，从而截留油边和眼球间的眼泪，并且睑脂使闭着的眼皮密闭。并且，它附着在眼泪表面，可以防止眼泪中的水过快的蒸发。功能障碍的睑板腺会导致眼睛干燥(由于没有油脂，眼睛中的水蒸发的过快)或睑炎；和其他的医疗状况。

已知，捕捉眼表的IR(infra-red，红外)图像来分析睑板腺的形态结构。对于健康的眼睛，腺体在空间宽度、面内伸长、长度等方面具有相似特征。另一方面，对于不健康的眼睛，成像的腺体显示不规则。因此，检测每个腺体和提取例如方向、宽度、长度、曲率等特征非常重要，以用于自动进行干眼症的诊断和风险评估。

图1(a)-(c)为眼表的三个IR图像样本。专家将图1(a)人工分级为“健康”，图1(b)为“中间”，图1(c)为“不健康”。在这个文件中，这样的图像分别称为“健康图像”、“中间图像”和“不健康图像”。IR图像具有一些使自动检测腺体区具有挑战性的特征：

腺体区和非腺体区的对比性低；

由光滑和潮湿表面引起的镜面反射；

由热成像引起的区域非均匀灰度分布；

眼表图像区的不规则。

即使成像条件存在多样性，腺体区与非腺体区相比具有更高的反射性。因此，属于腺体的图像区与相邻的非腺体区相比更亮。然而，由于以上提到的成像条件，传统方法例如局部阈值法不适于将图像划分为腺体区和非腺体区。

发明内容

本发明的实施例的主要目的在于，提供一种对包括多个睑板腺的眼部图像的自动处理，以识别睑板腺的位置和/或获得表征腺体的数值数据。这些数值数据可以用于将腺体分级。

为达到上述目的，本发明的实施例采用如下技术方案：

本发明的第一方面提出使用包含睑板腺区的眼部图像来获得指示睑板腺健康的等级，通过在眼部图像中使用来获得表征眼部图像显示的睑板腺的一个或多个数值参数；和使用一个或多个数值参数来自动确定等级。

等级可以用于筛选病人，以识别需要进一步详细检查的病人。它也可以用于提出要对病人实施的治疗方法。

所述数值参数最好包括以下至少一项：

(i)至少一个参数，所述参数表征由尺度空间变换获得的图像重要特征的尺度参数间的差异；和/或

(ii)至少一个参数，所述参数通过测量眼部图像中识别的、代表各自腺体的线来获得。所述测量可对单个线进行(例如线的长度)或涉及相邻线对(例如相邻线之间的距离)。

概括来说，本发明的第二方面提出使用伽柏(Gabor)滤波器作为局部滤波技术在眼部图像中识别睑板腺。

伽柏滤波器在形状、局部空间支持和取向上的参数化尤其适于检测睑板腺。

附图说明

现将参考下面附图，仅仅出于示例的目的，描述本发明的实施例，附图中：

图1由图1(a)-1(c)组成，示出三种捕获的眼表的IR图像；

图2由图2(a)-2(f)组成，示出六个伽柏函数的表示形式；

图3由图3(a)-3(f)组成，包括图3(a)，图3(a)为图1(a)的一部分，包括图3(b)-(f)，图3(b)-(f)阐明使用本发明实施例处理图3(a)图像的各个阶段。

图4由图4(a)-4(f)组成，包括图4(a)，图4(a)为图1(a)的一部分，包括图4(b)-(f)，图4(b)-(f)阐明使用本发明实施例处理图4(a)图像的各个阶段。

图5由图5(a)-5(f)组成，包括图5(a)，图5(a)为图1(a)的一部分，包括图5(b)-(f)，图5(b)-(f)示出使用本发明实施例处理图5(a)图像的各个阶段。

图6由图6(a)-6(f)组成，包括图6(a)，图6(a)为图1(a)的一部分，包括图6(b)-(f)，图6(b)-(f)阐明使用本发明实施例处理图6(a)图像的各个阶段。

图7由图7(a)和图7(b)组成，分别示出在进一步实施例中直方图均衡化前和后的眼部图像；

图8由图8(a)和图8(b)组成，图8(a)示出健康图像的尺度不变特征变换(尺度空间)点，图8(b)为不健康图像中的尺度空间点；

图9为健康和不健康图像的群体的尺度不变性分布和香农熵；

图10示出本发明进一步实施例中线特征的提取；

图11示出本发明进一步实施例使用的用于从像素集群中提取连续的线的方法；

图12示出由图11的方法提取的健康和不健康图像中线的总长度和线的数量的分布；

图13为健康和不健康图像中线的长度的散点图和长度的标准差；

图14由图14(a)和14(b)组成，示出叠加在眼部图像的线，包括表示腺体的线和虚假的线；

图15示出线端离图像边缘最大和最小距离的分布；

图16为图21的变型图，还包括中间图像的点；

图17为本发明第一实施例提供的流程图；

图18为本发明进一步实施例提供的流程图。

具体实施方式

第一实施例

1.1伽柏函数

本发明实施例是一种检测睑板腺的方法，其利用了2D(two-dimensional，二维)伽柏函数族。已知使用伽柏函数作为细胞的感受野函数，模拟简单细胞的空间总和属性[1]。使用伽柏函数的改进参数化将实验数据中发现的限制考虑在内[2、3]。假如一个二维视觉域Ω上的一点(x,y)(即)上存在一个光脉冲。伽柏函数由G_λ,θ,ψ(x,y)表示，其为一个实值数(real valued number)(即G_λ,θ,ψ(x,y)∈R)。伽柏函数由下式表示[2]：

$G_{\mathrm{\λ},\mathrm{\θ},\mathrm{\ψ}}(x,y)=\exp (-\frac{{\tilde{x}}^2+{\mathrm{\γ}}^2{\tilde{y}}^2}{2{\mathrm{\σ}}^2})\cos (2\mathrm{\π}\frac{\tilde{x}}{\mathrm{\λ}}+\mathrm{\ψ})-G_{\mathrm{\λ},\mathrm{\θ},\mathrm{\ψ}}^{\mathrm{DC}}$

其中

$\tilde{x}=-(x-x_0)\cos (\mathrm{\θ}-\mathrm{\π}/2)+(y-y_0)\sin (\mathrm{\θ}-\mathrm{\π}/2),$

$\tilde{y}=-(x-x_0)\sin (\mathrm{\θ}-\mathrm{\π}/2)+(y-y_0)\cos (\mathrm{\θ}-\mathrm{\π}/2).$

且为由余弦函数引起的直流项，直流项

$G_{\mathrm{\λ},\mathrm{\θ},\mathrm{\ψ}}^{\mathrm{DC}}=\underset{\mathrm{\Ω}}{\∫\∫}{G}_{\mathrm{\λ},\mathrm{\θ},\mathrm{\ψ}}(x,y)\mathrm{dxdy}$

为G_λ,θ,ψ减去偏压。

不失一般性的，在实施例中假设伽柏函数集中在感受野的坐标平面上。因此，x₀和y₀不能用于索引(index)感受野函数。参数σ,γ,λ,θ和ψ在下面解释。

感受野的大小取决于高斯因子的标准差σ。参数γ的范围为0.23-0.92(即γ∈(0.23,0.92))[2]，被称为空间纵横比。它决定感受野的椭圆率。在下面的实验结果中使用值γ＝0.5，由于此值为常数，参数γ不能用于索引感受野函数。参数λ为波长，1/λ为余弦因子的空间频率。比值σ/λ决定空间频率带宽，并因此决定在感受野中观察到的如图2所示的平行兴奋和抑制条纹区的数量(以下给出解释)。根据公式(1)，具有脉冲响应的线性滤波器的半响应空间频率带宽b∈[0.5,2.5](八度(octave)内)[2]为以下的比值σ/λ的函数[2]:

$b={\log }_2\left(\frac{\frac{\mathrm{\σ}}{\mathrm{\λ}}+\frac{1}{\mathrm{\π}}\sqrt{\frac{1n\left(2\right)}{2}}}{\frac{\mathrm{\σ}}{\mathrm{\λ}}-\frac{1}{\mathrm{\π}}\sqrt{\frac{1n\left(2\right)}{2}}}\right)---\left(3\right)$

或者相反的

$\frac{\mathrm{\σ}}{\mathrm{\λ}}=\frac{1}{\mathrm{\π}}\sqrt{\frac{1n\left(2\right)}{2}}\frac{2^b+1}{2^b+1}---\left(4\right)$

在实施例中值b＝1.0，由于这个值为常数，当给定λ时，参数σ可以由公式(4)求出，因而不能用于索引感受野函数。角参数θ∈[0,π)沿逆时针方向决定从x轴的优选取向。参数ψ∈(-π,π]为相位偏移，其决定G_λ,θ,ψ(x,y)关于原点的对称性：当ψ＝0和ψ＝π时，它是对称的(或偶数)，当ψ＝-π/2和ψ＝π/2时，它是反对称的(或奇数)；所有其他情况为不对称的混合。

图2为模拟一简单细胞的感受野剖面的伽柏函数的强度图，伽柏函数的参数为θ＝π/4、ψ＝0.0及在256×256像素网格上(a)λ＝10像素；(b)λ＝20像素；(c)λ＝30像素；(d)λ＝40像素；(e)λ＝50像素；(f)λ＝60像素。比背景更亮和更暗的灰色阶分别指示函数取正值和负值的区域。亮椭圆表明(典型的)感受野的边界。在该边界外，函数取极小值。

伽柏函数取以上参数值时，我们可以计算对于输入2D图像I的响应I_λ,θ,ψ为

I_λ,θ,ψ＝I*G_λ,θ,ψ  (5)

其中，*表示二维卷积，使用傅里叶变换(F)，即I_λ,θ,ψ＝F^-1(F(I)F(G_λ,θ,ψ))，可以有效的计算出公式(5)，其中F^-1为反傅里叶变换。

1.2使用伽柏滤波器提取的特征

图2中显示的伽柏函数的实现可被用于模拟被非腺体区包围的腺体的局部结构。也就是说，中间的主瓣代表腺体，主瓣两边的旁瓣代表非腺体区。参数λ可以用于评估空间宽度，参数θ可以用于次级腺体结构的局部取向。使用ψ＝0.0。不失一般性的，参数ψ不用于索引，除非另有说明。

参数λ从一个有限集{λ_i}中取离散的整数值，并可以通过连续的腺体区和非腺体区的预期空间宽度来估算。同时，可以预计次级腺体结构可以有[0,π)之间的任意取向。然而，不可能测试出所有可能的取向，参数θ通过下式进行离散：

${\mathrm{\θ}}_i=\frac{(i-1)\mathrm{\π}}{\mathrm{N\θ}},i=\mathrm{1,2},...,N_{\mathrm{\θ}},---\left(6\right)$

其中N_θ为离散取向的总数。

对于每个像素的修正值λ和θ的粗略评估，伽柏滤波器在腺体区的响应为正，在非腺体区为负。为了证明以上观点，使用如图3所示的来自于图1的IR图像的次级区域。下面我们将阐明使用图3(a)、数值λ＝40像素、N＝180获得图3(f)示出的二值化(binarized)滤波器响应的步骤。图3(b)为对比增强IR图像，其中图3(f)的区域的边界像素被重叠，四个像素位置分别标为A、B、C和D。像素A和B位于腺体区，像素C位于非腺体区。像素D位于腺体区和非腺体区的接界。

图3(c)为图3(b)中的四个像素的伽柏滤波器响应。图3(c)相对于变化的θ作图，图3(c)记录了落入腺体和非腺体区的不同像素位置的伽柏滤波器的响应。当伽柏滤波器的响应信号为+1时，实现了像素A和B的伽柏滤波器响应最大绝对值，与此同时，像素C的伽柏滤波器响应最大绝对值为-1。因此，对于被非腺体区包围的腺体区的像素，正最大值和负最小值的量级间有很大的不同。类似的，对于被腺体区包围的非腺体区的像素，负最小值和正最大值之间的量级间有很大的不同。然而，当像素位于非常靠近腺体和非腺体之间的区域时，例如图3(b)中的像素D，正最大值和负最小值的量级间的差距很小，如图3(c)所示。并且，清楚的是，在腺体区伽柏响应的平均值标定为1，非腺体区为-1。

当给定λ时，平均伽柏响应用下式计算：

${\hat{I}}_{\mathrm{\λ}}=\frac{\underset{0}{\overset{\mathrm{\π}}{\∫}}{I}_{\mathrm{\λ},\mathrm{\θ}}\mathrm{d\θ}}{\underset{0}{\overset{\mathrm{\π}}{\∫}}\mathrm{d\θ}}=\frac{1}{\mathrm{\θ}}\underset{0}{\overset{\mathrm{\π}}{\∫}}{I}_{\mathrm{\λ},\mathrm{\θ}}\mathrm{d\θ},---\left(7\right)$

可以近似为

${\hat{I}}_{\mathrm{\λ}}=\frac{1}{N_{\mathrm{\θ}}}\underset{i=1}{\overset{N_{\mathrm{\θ}}}{\mathrm{\Σ}}}{I}_{\mathrm{\λ},{\mathrm{\θ}}_i}.---\left(8\right)$

图3(d)示出了当λ＝40像素且N_θ＝180时，由公式(8)计算的图3(a)的输入图像的平均伽柏滤波器响应图3(e)给出了相应的表面图。从图3(d)和图3(e)中可以明显得出在腺体区，滤波器给定高的正响应，在非腺体区，滤波器产生低的负响应。因此，使用此观察结果，我们可以使用滤波器响应的标志很容易的从非腺体区中分割出腺体区，例如

${\hat{B}}_{\mathrm{\λ}}(x,y)=H\left({\hat{I}}_{\mathrm{\λ}}(x,y)\right).---\left(9\right)$

其中H(a)∈{0,1}为香农函数，其定义为

图3(f)示出了根据公式9的二值化伽柏滤波器响应总体的，分割结果看起来令人满意。然而，可以观察到，像素D处的两个腺体区被合并在一起，其中两个腺体间的间隔在图3(e)所示的表面图中是明显的。这可以由于以下两种原因引起：1)λ值对于解决信号分离而言太大；和/或2)像素D落入不确定的区域，其中，没有足够的信息可以分离两个区域或者此区域属于眼表的其他部分。在后面一种情况，需要进一步的分析。然而，前种情况可以通过分析不同A值处的图像来解决。例如在图4中，使用A＝20像素和N_θ＝180来分析图3(a)中的相同输入图像。图4(a)与图3(a)相同，图4中每一其他部分分别对应于图3中的各自的图像。清楚的是，像素D正确分割的代价为其他区域的不正确分割。

图5示出了λ＝30像素和N_θ＝180时的相应结果，其中每个图像分别对应于图3的图像。

利用不同的λ值获得的结果表明在空间细节保留和降噪之间产生不同的权衡。尤其，具有较低λ值的图像更容易受到噪声干扰，同时保留图像内容的更多细节。相反的，具有较高λ值的图像不容易受到噪声的干扰，但是会牺牲较多的图像降级细节。

对于λ∈{λ_i}中的每个不同值的平均伽柏滤波器响应都用于输入图像的每个像素(x,y)的矢量表示f_x,y中，为

其中，N_λ为集{λ_i}中的点，i为正整数(即i∈Ζ⁺)，且

$f_{x,y}\left(i\right)=\frac{I_{{\mathrm{\λ}}_i}(x,y)}{\sqrt{\underset{k=1}{\overset{N_{\mathrm{\λ}}}{\mathrm{\Σ}}}{I}_{{\mathrm{\λ}}_k}^2(x,y)}}---\left(12\right)$

公式(12)中的分母用于补偿由图像的不同部分的照明差异引起的波动。图6示出了N_λ＝16和N_θ＝180时由公式(12)得出的伽柏滤波器响应。图6(a)与图3(a)、4(a)和5(a)相同。下面我们将会解释使用图6(a)获得的图6(f)示出的局部滤波器响应的步骤。图6(b)为对比增强IR图像，其中图6(f)的区域的边界像素被重叠，且四个被标为A、B、C、D的像素位置与图3-5中的一样。图6(c)示出了图6(b)中四个像素通过公式(12)得到的伽柏滤波器响应。图6(c)绘出了落入腺体区和非腺体区的不同的像素位置的伽柏滤波器响应。从图6(c)所示的特征向量可知，在腺体区的特征向量是正的，非腺体区是负的。然而，它们在像素D的正值和负值间波动，平均特征由下式计算

$\hat{F(x,y)=\frac{1}{N_{\mathrm{\λ}}}}\underset{i=1}{\overset{N_{\mathrm{\λ}}}{\mathrm{\Σ}}}{f}_{x,y}\left(i\right)---\left(13\right)$

如图6(d)和(e)所示，图中腺体区和非腺体区的区别很明显。使用公式(13)最后的二进制图像由下式计算

$\hat{B}(x,y)=H\left(\widehat{F(x,y)}\right)---\left(14\right)$

其中，H(·)由公式(10)定义。图6示出由公式(14)得到的二进制化结果，且边界被叠加到图6(b)中的输入图像上。清楚的是，λ的不同值的合并信息提高腺体区和非腺体区的区别。

因此，图17为第一实施例的步骤。在步骤1中，λ值被初始化(即i设置为第一个值，以便设置λ_i)。在步骤2中，θ被初始化。在步骤3中，λ值和θ用于执行伽柏滤波器变换。针对每个可能θ值重复该步骤，得到的结果用于步骤4以得到对应每个λ值的I_λ值。该结果针对可能的λ值被迭代，得到的结果使用公式(13)和(14)获得步骤5中的

第二实施例

本发明的第二实施例主要目的是提供一种分级受试者(subject)的方法，即将他分为至少两种种类中的一种，例如“健康”、“不健康”或“中间”。

图18描述了第二实施例的所有方法。单个眼部图像用于获得指示一个或多个表征图像是健康的或不健康的数值参数(“特征”)。请注意，并不是下面描述的所有的数值参数都在实现实施例的过程被收集，但是收集的参数最好超过一个，在这种情况下，数值参数由自适应学习系统(如支持向量机(SVM)，其一直受到监督学习)结合，来产生指示图像是否健康的输出。

2.1使用直方图均衡化增强对比度

原始图像的对比度差，为了提升对比度，我们使用一个称为直方图均衡化的标准技术。图7(a)示出原始图，作为对比，图7(b)示出提高对比度后的图像。在经过直方图均衡化的图像上执行第二实施例的操作。

2.2尺度空间香农熵特征

2.2.1尺度空间变换

第二实施例利用尺度空间香农熵特性使健康图像区分于不健康图像。这个概念来自于[4]中描述的著名尺度不变特征变换(SIFT)方法。简而言之，在图19中的步骤20，实施例在图像上定位关键点，称为尺度空间点。每个尺度空间点由一个具有3个元素的矢量(x,y,s)表示(注意，与之相比，SIFT变换使用了进一步的129个元素，其并没有被本实施例采用)。x和y为图像上的尺度空间点的笛卡儿坐标，s称为它的尺度。通过下面的“尺度空间变换”找到尺度空间点。

-在不同距离尺度s处对眼部图像(x,y)与高斯过滤器卷积；

-寻找具有不同s的后续图像对之间的差异的最大值(x,y,s)，以形成候选关键点；

-去除与下一个邻近点相比低于一个阈值的候选关键点。

图8示出在两个图像即健康图像(图8(a))和不健康图像(图8(b))上，尺度空间点看起来的样子。每个尺度空间点被代表为一个圆，圆的水平横条(即半径)代表它的尺度。

实施例利用如图8所示的局部区域范围内的观察(由方框显示)，健康图像中的圆圈为相似的尺寸，而不健康图像的圆圈为不同的尺寸。这是因为在健康图像中，存在有相似厚度的分布均匀间隔带(即腺体)，尺度空间变换获得了这种图样。另一方面，不健康的图像没有这种模式(即没有腺体)，圆圈的尺寸是非常不同的量级。

在图18中的步骤21，实施例产生一种圆圈的尺度差异的数值测量，且实施例利用健康图像中尺度局部分布均匀而不健康的图像不均匀这样的事实来区分两种等级。可以测量此种均匀性的数学函数为著名的香农熵，现在我们将对其进行讨论。

2.2.2香农熵

香农熵的定义为

$S=-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i\ln p_i---\left(15\right)$

其中p_i为事件i的概率。p_i必须为归一化的，即

我们使用的香农熵的一个重要性质为有且只有p_i为均匀分布时，即p_i＝1/n，香农熵为最大。其证明参考[5]。

SIFI点的尺寸可以以下方式与概率分布建立关联。首先，选择一个尺度空间点，考虑它的n个最近的邻近空间点。我们定义关于“中心”尺度空间点的第i个最近的邻近概率为

$p_i=\frac{{\left(s_i\right)}^2}{{\mathrm{\Σ}}_{i=1}^n{\left(s_i\right)}^2}---\left(16\right)$

其中i＝1,...,n标记n个最邻近空间尺度点，s_i为第i个空间尺度点的尺寸。

公式(16)中的分母保证0<p_i<1，且其分布为归一化的。p_i的含义为第i个邻近点的圆圈面积与所有邻近点的总面积的比值。

参考图8，我们看到对于图8(a)中的方框，所有的圆圈都近似为相同的尺寸，所以p_i也近似相同，因此在蓝框内计算的熵应该接近于最大值。另一方面，对于图8(b)，由于圆圈不均匀，由公式(16)计算得到的p_i也不均匀，因此在蓝框内计算的熵与(a)相比将会低。

2.2.3计算图像的香农熵的算法

为了比较两个不同图像的熵，我们需要计算整个图像的香农熵。以下为其算法。

1、获得一个图像的所有尺度空间点(使用标准技术)。假设尺度空间点的总数为M，用α表示一个空间尺度点，其中α＝1,...,M。

2、对于一个α，定义它的n个最近邻近点(通常n＝20)。(假如其中n个最近邻近点没有被明确定义(例如，因为有3个相距距离完全相同的尺度空间点和19个相近的尺度空间点；或者更普遍的说，如果对于最小距离d，在α中至少有n个点不远于d，α中不远于d的空间尺度点的个数为m，因此m其大于n)，那么算法可以随机的选择这些m个点中的n个，或者可替代的，使用这些点中的全部m个点。)

3、如果{s₁ ^α,...,s_n ^α}为这些n个最邻近点的尺寸，令s^α＝∑_i＝1 ⁿ(s_i ^α)²，且令p_i ^α＝(s_i ^α)²/s^α。

4、计算α的熵：S^α＝-∑_i＝1 ⁿp_i ^αlnp_i ^α。

5、对于所有的空间点即α＝1,...,M，重复步骤1到步骤3。

6、整个图像的香农熵为平均值

$\stackrel{\&OverBar;}{S}=\frac{1}{M}\underset{\mathrm{\&alpha;}=1}{\overset{M}{\mathrm{\&Sigma;}}}{S}^{\mathrm{\&alpha;}}---\left(17\right)$

2.2.4结果

为了查看尺度空间香农熵的特征是否可以区分健康的和不健康的睑板腺成像图像，我们将一些健康和不健康的图像人工分级，计算其香农熵S。结果在图9显示。每个图像用一个点表示。横轴表示熵S。纵轴只包括可视化的情况。健康图像绘制为更亮的点，不健康的图像为各自更暗的点。结果表明，沿着香农熵维度，健康集群可以与不健康集群很容易的区分。

2.3线特征

实施例使用进一步的方法从图像中提取其他特征，我们称之为线特征。如所说明的，健康图像最显著的特征为垂直腺图样，如图10中左边最高板块所示。图10的底部两个板块为本实施方式获得的指示这些腺体的像素的集群。其过程可以概括如下。

1、提取沿着明亮区和黑暗区(即腺体区域)的像素。这是图18的步骤22。

2、将像素分组为类似于线图样的“原始”集群。这是图18中的步骤24。

3、对每个集群进行形态运算，形成集群中的一条线。这是图18中的步骤24。

4、获得这些线的数值特性(例它们的长度和/或曲率)作为分类的特征。这是图18中的步骤25。

接下来，我们将详细描述上面提到的每个步骤。

2.3.1提取局部最小和最大像素

为了提取局部最小值，我们首先使整个图像平滑(使用著名的程序函数的OpenCV(开放源代码计算机视觉类)库的中的cvSmooth(平滑)算法)以使强度值可以平滑的从一个像素变化到另一个。接着，将此强度做成一种“势能面”，我们找到面的最小值。

这样做的一种方式为执行梯度下降以达到面的最小值。具体的步骤如下所示。

1、考虑图像的水平截面，即一行。定义第i个像素的强度为I[i]。

2、考虑位于第i个像素的一个“点”。

3、计算点上的“力”为

f[i]＝-(I[i+1]-I[i-1])  (18)

(其仅仅是梯度公式的简化版。)如果f[i]>0，将点移到i+1；如果f[i]<0，将点移动到i-1，如果f[i]＝0，点停留在i。注意，由于强度通常引用整数值(例范围0到255的整数值)，我们可以预测最小值处f[i]＝0。

4、重复步骤3，直到点停止移动，即f＝0。定义点的位置为i^*。

5、重复步骤2和3，使粒子从行的每个像素开始，获得对应于每个开始像素i的i^*。所有的i^*设置为行的最小值。

6、继续下一行，重复步骤2-5。

为了获得最大值，我们简单的重复以上步骤，但是对于像素强度使用-/代替/。

另一种替代算法是沿着行扫描，制作一个满足下列条件的那些点/像素i的列表:对于最大值：I[i-1]<I[i]>I[i+1]；对于最小值：I[i-1]>I[i]<I[i+1]。

2.3.2聚集像素

获得所有最小点和最大点后，我们需要聚集它们，使得，对于第一次逼近而言，每个集群都有线类似的形状。这将促进集群转换成连续的线这样的下一步骤。为了完成我们的集群，我们使用称为先进先出(FIFO)的著名算法。具体的，我们使用10像素的阈值，这意味着只有集群中一些其他像素中少于10个像素的像素才能被组合在一起。这个过程之后，最小值(最大值)被组合，如图10的两个底端框图所示。接着使用下一节描述的方法将每个集群被转换为一条连续的曲线。

2.3.3将集群转换为曲线

在这个阶段，每个集群类似于一条线，但它仍然没有用，因为它可能被打破，且它也不是一个像素厚(pixel thick)。在这里，我们基于众所周知的方法提出一种算法将集群转换为曲线。虽然实施方式中使用的每个方法都为无关的方法，本实施例结合这些方法将像素集群转换为连续的曲线。

图11为图18中的步骤24的第二个步骤，其描述如下所示。

1、选择一个集群，将它放入一个图像中。具体来说，我们设置所有属于那个集群的像素作为前景(图1中的黑色)，所有其他像素作为背景(即白色)。如图11中的第一个框图所示。

2、接着，我们运用cvDilate(膨胀函数)(OpenCV中的另一个算法)由一个像素来加厚集群。其目的是将所有像素合并到一个连接片。应用一次cvDilate后，我们检查看集群目前是否由单个的连接组件组成。如果是，我们进入到下个步骤，否则，我们再次应用cvDilate直到获得一个单一连接组件。

3、前一步骤可能产生一个包含背景的“岛屿”的连接组件，由图11中的第二个板块中的圆圈标出。必须消除这些岛屿因为它们会在最后步骤中，在连续的线上产生“线圈”。为了消除它们，首先我们使用cvFloodFill(漫水填充)(OpenCV中的另一个算法)来填充背景，将这些岛屿的位置暴露为剩余的白色像素。接着，我们回到使用cvFloodFill前的图像，设置这些岛屿为前景(即黑色)。这会产生一个不包含这些问题岛屿的连接组件。

4、接着，我们使用一个标准细化算法(thinning algorithm)[6]来细化(thin)组件直到得到一个像素厚。

5、连接组件现在成为一个具有很多枝的“树”。实施例减掉了所有的旁枝。为此，我们首先定位和计算树上端点(即结束点)的数量。然后，我们使用一个标准的剪枝算法(pruning algorithm)[7]用修剪因子来剪掉所有旁枝，即长1像素的枝被删除。接着，我们再次计算端点的数量，如果不是2，我们对修剪因子加1，然后再次修剪。重复该过程直到只剩下两个端点，这意味着树已经没有剩余的枝了。

6、这样就完成了将一个集群转换成一个连续的线的过程。接着，对于另一个集群，我们重复步骤1-5。

2.3.4健康和不健康图像的特征选择和分类

现在，我们展现在图18的步骤25中获得的线推导出来的特征(即数值参数)，来帮助我们区分健康和非健康图像。特征为:

1、线的数量。

2、线的总长度。

下面所有的变换特征可能会被使用。

3、势能。

4、左右距离。

5、扭曲。

线的数量和总长度对于前两个特征，我们选择“线的数量”和“线的总长度”。这个想法的解释在图12中。对于健康图像，只有很少的线，且所有的线都很长；对于不健康图像，存在更多的线，且线都很短。如图12的底部框图所示，如果我们利用这两个特征，可以区分健康和不健康图像。

势能下一个特征称为“势能”。首先，我们找到最小值(最大值)点的最邻近点。接着我们计算它的势能。除了距离约为50像素时(约为条间距)势能为负值，其他距离下势能均为0。因此，对于健康图像，我们预计图像的总势能为负值，而对于非健康的图像，我们预计它近似于0。

左右距离下一个特征称为左右距离。我们从一个点开始，沿垂直于该点切线的两个方向移动。对于每个方向，我们计算该实施例首次遇到另一个点(图中的d₁和d₂)的距离。点在d₁-d₂平面的分布可以用于构造一个直方图。这个直方图的前两个组件为分离健康和非健康图像的很好特征。

扭曲我们最后一个特征称为扭曲，它是基于不健康图像的线比健康图像的线扭曲的观察。为了量化扭曲度，我们将它定义为线的端点间的距离(“直线长度”)和线的长度(“曲线长度”)的比值。对于每个图像，可以计算直方图用于其线的扭曲的分布。接着，我们进一步粗糙颗粒化直方图使其成为两个箱子(bin)，使其作为特征空间。

2.4第二实施例中的变化

以上讨论的特征并不是利用步骤25从步骤24获得的线中推导出的唯一特征。另两个更进一步的可能特征为图像中的所有线的长度的平均值，和线的长度标准偏差。图13示出基于这两个特征的图像的散点图。每个点代表一个图像。形状为圆圈状和钻石状的点分别代表健康和不健康图像。清楚的是，使用如支持向量机的标准技术，图像的种类可分为两种。这种方法已经在我们最近发表的生物医学光学期刊2012年17卷第8期0860008(Journal of Biomedical Optics17(8)0860008,(2012))中报道。

此外，在用于获得线的技术中许多变化是可能的，也可以做出改进。虽然上面给出的从图像中提取线特征的算法对健康和不健康图像的极端情况的分类是有效的，但是对于评估中间阶段可能不太有效，且它对于降低噪声有一定的价值。这里的噪声指由图18的算法提取的但不是睑板腺导致的虚假线。这些线中的大部分都是由于眼睫毛或图像边缘的阴影导致的。这些被排出，以便消除我们从线中计算得到的后期统计数据中的噪音。图14(a)中显示了一个例子。由图18的方法推导的线被标记。标记为X的线是由睑板腺引起的，但是接近边缘的标记为Y的线为虚假的，是由于像素强度的不均匀引起的，其最好是在实施步骤25前排除，即从基于线的数据的后续计算中排除。随后，我们描述用于自动检测虚假线Y的算法。

虚假线Y有两个重要的特性。它们比大部分线都短且它们接近于边缘。然而，如图14(b)所示，将所有短的线都排除是不合适的，因为它们是不健康和中间图像的特征。我们看到，由于位于中心的腺体模式的破裂，在图像的中心有许多短虚线。这些必须被保留。

排除这些假线有两个阶段。在第一阶段中，我们用以下方式处理每个单独的线。

(1)标定图像的顶端、底端、左端和右端(例如分别为1、2、3、4)。对于线的每个像素，检查像素接近哪端，确定像素到最近端的距离。称这个号为d。

(2)对于每条线，得到具有最小d和最大d的像素，分别称它们为d_min和d_max。

(3)计算图中所有线的d_min和d_max。在一个二维空间中将每条线作为一个点画出来，其中x轴为d_max-d_min，y轴为d_max。其结果的散点图如图15所示。

所有的线都落入二维图的三个之一分区。代表腺体的线会落在顶端右手区。这是因为它们通常包含一个接近边缘的像素，因此有一个小的d_min，但是也包含一个接近中心的像素，因此有一个大的d_max。结果，与腺体相关的线将会有大的d_max-d_min和d_max。

第二种线属于落入顶端左手区的线。它们是破裂腺体线的残余线。它们很短且接近于中心。由于它们接近于图像的中心，所以它们的d_max是大的。但是因为它们很短且接近于中心，具有d_min的像素通常也靠近于具有d_max的像素，且d_min近似等于d_max，给出一个小的d_max-d_min。这解释了为什么它们位于顶端的左手角。当实施步骤25时，这些线也应该被包括在内。

这些虚假线Y为位于接近散点图的起源的线。它们是短的，因而沿着线的每个像素将会接近于边缘，因此d_max将会很小，d_max-d_min也会很小。

为了选择虚假线，我们使用一个阈值方法。符合d_max-d_min<50和d_max<50的任意线都会被定义为虚假线，因此其在步骤25实施前被排除。

2.5工业实用性

图13表明使用图像中所有线的长度的平均值和线的标准偏差可以很好的区分健康和不健康图像。这个结果和以上给出的其他的实验结果表明，第二实施例可以作为一个成功的工具用于将受试者识别为这些类别，甚至使用单一的眼部图像就可以实现。被发现具有不健康的眼睛的受试者可以接受进一步检查，或治疗。众所周知，用热毛巾、眼罩或者叫做睑板蒸汽的专门加热装置的温暖压缩来治疗睑板腺功能障碍患者。此外，抗炎药物如强力霉素、阿奇霉素和环孢菌素可能对其有帮助。另外，患者可以开始外敷抗生素类固醇药膏。除了计数润滑油，尤其是含有可以补充的眼泪脂质的那些脂质，因为具有这些情况的病人可能存在泪膜异常脂质层。这种情况下的最新进展包括具有或不具有类固醇的腺内注射睑板腺的机械探测，以及称为热力脉冲振动仪器的办公室睑板腺脉搏治疗。由于这些治疗中的很多可能倾向于依从性差，或者涉及并发症或者昂贵，因此拥有诊断和评估睑板腺功能障碍的客观方法是有帮助的。

图16使用与图13相同的两个特征，但是也包括中间图像的点，其表明对于这两个特征，中间等级与健康和不健康的等级很明显的重叠。因此，使用其他特征例如以上描述的那些来帮助区分中间图像是可取的。例如，意图使用上面描述的用于区分健康和不健康图像的分类器来首次给出一个初步的分类，将图像分成两类(i)健康/中间图像或(ii)不健康/中间图像。区分健康/中间图像的特定特征可以接着用于进一步在健康/中间类别中将健康与中间类别分离。对于不健康/中间类别，方法相似。

参考文献

[1]J.G.Daugman,"Uncertainty relation for resolution in space,spatialfrequency,and orientation optimized by two-dimensional visual corticalfilters,"J.Opt.Soc.Am.A,vol.2,no.7,pp.160-1169,Jul1985.

[2]N.Petkov and P.Kruizinga,"Computational models of visualneurons specialised in the detection of periodic and aperiodic oriented visualstimuli:bar and grating cells,"Biological Cybernetics,vol.76,pp.83-96,1997.

[3]P.Kruizinga and N.Petkov,"Nonlinear operator for orientedtexture,"IEEE Transactions on Image Processing,vol.8,no.10,pp.1395-1407,Oct1999.

[4]D.G.Lowe,"Object Recognition from Local Scale-InvariantFeatures,"The Proceedings of the Seventh IEEE International Conference onComputer Vision,vol.2,pp.1150-1157,1999.

[5]A.I.Khinchin,Mathematical Foundations of Information Theory.Dover Publications,1957.

[6]L.Lam,S.Lee,and C.Suen,"Thinning methodologies-acomprehensive survey,"IEEE Transactions on Pattern Analysis and MachineIntelligence,vol.14,pp.869-885,1992.

[7]A.Niemisto,V.Dunmire,O.Yli-Harja,W.Zhang,and I.Shmulevich,"Robust quantification of in vitro angiogenesis through imageanalysis."IEEE Transactions on Medical Imaging,vol.24,no.4,pp.549-553,2005.

[8]K.K.Nichols,G.N.Foulks,A.J.Bron,B.J.Glasgow,M.Dogru,K.Tsubota,M.Lemp,and D.A.Sullivan,The International Workshop onMeibomian Gland Dysfunction:Executive Summary,Ophthalmol.Vis.Sci.March30,2011vol.52no.41922-1929.

Claims (15)

1.一种由计算机设备执行的、使用包含睑板腺的区域的眼部图像来获得指示睑板腺的健康等级的方法，所述方法包括：

(i)自动获得表征所述眼部图像中显示的睑板腺的一个或多个数值参数；

(ii)使用多个数值参数中的一个来自动确定等级。

2.根据权利要求1所述的方法，其中，所述获取所述一个或多个数值参数的步骤包括：

在所述图像中产生一个关键点，每个关键点与各自的距离尺度值相关；

获得所述数值参数中的一个，作为指示所述关键点的所述尺度值的差异的参数。

3.根据权利要求2所述的方法，其中，所述的数值参数被计算为关键点的各自数值S的的平均值，每个关键点的数值S的计算基于最接近所述关键点的其他关键点的子集。

4.根据权利要求3所述的方法，其中，所述数值S根据以下表达式计算

$S=-\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{p}_i\ln p_i$

其中，i＝1,...,n表示子集的n个关键点，p_i由下式给出

$p_i=\frac{{\left(s_i\right)}^2}{{\mathrm{\&Sigma;}}_{i=1}^n{\left(s_i\right)}^2}$

其中，s_i为n个关键点的尺度值。

5.根据前述任一项权利要求所述的方法，其中，通过以下步骤推导出至少一个所述数值特征：

获取眼部图像中的线，其指示各自的睑板腺；

获取至少一个由使用线而测量到的数值特征。

6.根据权利要求5所述的方法，其中，所述获取图像中的线的步骤包括：

提取沿图像的明亮区和黑暗区分布的像素；和

将像素分组为集群。

7.根据权利要求6所述的方法，其中，所述获取线的步骤还包括对每个集群进行形态运算以形成来自集群的一条线。

8.根据权利要求5-7中任一项所述的方法，包括线排除步骤：识别所述线中与腺体无关的线，将它们从考虑的线中除去。

9.根据权利要求8所述的方法，其中，所述线排除步骤的实现包括：确定每个线到图像边缘的最小距离d_min，到图像边缘的最大距离d_max，和使用所述距离识别所述线中与腺体不相关的线。

10.根据权利要求9所述的方法，其中，所述与腺体不相关的线被确定为那些d_max和d_max-d_min低于各自阈值的线。

11.根据权利要求5到10中任一项所述的方法，其中，使用线获得的数值参数包括以下中任意一个或多个：

(i)线的数量；

(ii)线的总长度；

(iii)线的平均长度；

(iv)线的长度的标准差；

(v)通过多条线得到的值，所述值通过加和每条线上的作为该条线上的点到另一条线上的点的距离的函数获得的值而得到；

(vi)通过多条线得到的值，所述值通过加和每条线上的通过测量从该条线沿其切线到到另一条所述线的至少一个距离获得的值而得到；

(vii)指示线的端点间的距离和相应线的长度的比值的值。

12.一种划分包括睑板腺区的眼部图像的方法，所述方法包括：

在图像的多个位置中的每一位置处：

(i)将以该位置为中心的图像进行伽柏函数变换，在图像中至少用一个尺度因子λ和一个方向θ来表征；

(ii)将各θ值的伽柏函数变换加和，来形成一个强度值和

(iii)使用执行阈值步骤，来得到一个代表相应位置是否对应于睑板腺位置的二进制值。

13.根据权利要求12所述的方法，其中，步骤(iii)中，值为λ的多个值的总和，且结果为一个阈值。

14.一种分析眼部图像的计算机设备，包括处理器，数据存储装置，所述数据存储装置存储计算机指令，所述计算机指令当被所述处理器执行时运行，使所述处理器执行根据权利要求1-13任一项所述的方法。

15.一种计算机程序产品，包括程序指令，所述程序指令当被处理器处理时运行，使所述处理器执行根据权利要求1-13任一项所述的方法。