1 引言

光学相干层析成像(OCT)是一种新兴的生物医学光学成像方式,可对被测样品进行在体无损成像。OCT以其高分辨率、非侵入和无辐射等优点,在眼科临床中得到广泛应用^[1-2]。

角膜厚度的改变是角膜类疾病的重要表现因素。角膜厚度的精确测量可应用于干眼症、角膜炎、青光眼等角膜类疾病的诊断,指导治疗方案、手术方式^[3-7]。OCT成像设备可以快速、非接触地获取角膜的高分辨率断层成像,为角膜厚度的精确测量提供可能。临床上对角膜图像的处理一般有3种做法:1)医生直接观察图像中角膜形状的细微变化,凭借经验分析图像,判断有无病变;2)医生在角膜图像中选取几个点,程序利用这几个点计算出角膜的生物结构参数;3)使用数字图像处理技术,根据角膜的结构特征,分割出角膜边界,对其进行曲线拟合,得出角膜生物数据^[5-6,8-16]。手动测量需要医生花费较多的时间分析角膜图像和计算角膜的生物结构参数,自动测量方式可在不超过1 s的时间内计算出角膜的厚度。相比于手动测量方式,自动测量效率更高,受主观因素影响小,可重复性好^[5]。

在实际应用中,OCT系统容易受测试环境和被测样品的影响,导致图像中存在多种噪声;大多数频域光学相干层析成像(FDOCT)角膜成像仪采用的是远心(平行)扫描模式,导致在角膜顶点处常出现中心饱和伪影。这些问题都给角膜厚度的自动测量带来了困难^[7-9]。

为了能够自动分割角膜边界并获取角膜厚度,研究者们提出了几种不同的方法来分割角膜层边界。Feizi等^[17-19]提出了一种快速活动轮廓和二阶多项式拟合相结合的自动角膜分割方法,该方法需要将初始轮廓放置在接近期望灰度特征(通常是边缘)的位置,该轮廓在外力和内力作用下产生形变,直到收敛到感兴趣区域。该方法的处理过程有两个关键点:1)如何放置这一轮廓;2)如何使这一轮廓变形,并收敛到所需特征。这两个关键点的选择影响图像边缘轮廓的提取,同时也容易受到图像中噪声的干扰。Li 等^[20]在2007年利用增强的智能剪刀(Eis)实现半自动的分割方法,用户交互的分割方法需要最少的用户输入,以及能量最少的样条;该方法中,第一阶段应用形态学算子来增强角膜的突出结构特征,第二阶段使用半自动分割算法分割角膜的上下边界,该方法虽然具有较强的抗干扰能力,但无法对存在中心饱和伪影的角膜图像进行很好的分割。

舒鹏等^[7,21-24]提出了一种基于边缘检测和随机抽样一致性算法(RANSAC)的中央角膜厚度自动测量方法;该方法采用边缘检测算子提取图像中存在的边缘,然后利用随机抽样一致性算法对图像轮廓边缘进行曲线拟合,曲线拟合的精度由随机抽样一致性算法的迭代次数决定,迭代次数越多,拟合精度越高。但对于存在较多噪声和中心饱和伪影的图像,拟合效果并不理想,存在较大误差。

本文提出了一种基于边界跟踪、轮廓定位和曲线拟合的角膜厚度自动测量方法。在利用自适应中值滤波和伪影去除方法对图像进行预处理的基础上,使用边界跟踪算法获取角膜图像的边缘轮廓,根据角膜边缘轮廓点拟合出角膜上下边缘,得到拟合曲线方程,进而计算出角膜厚度。该方法对临床角膜OCT图像中常见的伪影、噪声和低信噪比区域具有很强的稳健性,对存在角膜中心饱和伪影的图像也能够实现精确测量,具有抗干扰能力强和测量精度高的优点。

2 基本原理

所提出的角膜厚度自动测量方法流程如下:

1) 利用自适应中值滤波和伪影去除方法对原始角膜OCT图像进行预处理;

2) 对预处理后的图像进行二值化处理,采用边界跟踪算法获取角膜上下边缘轮廓点;

3) 根据获取的角膜上下边缘轮廓点,利用最小二乘法拟合出角膜的上下边缘,得到拟合曲线,求出上下边缘任意一点在整幅角膜图像中的坐标;

4) 将角膜上下边缘坐标点代入角膜厚度计算公式,计算出角膜厚度。

该方法不仅可以准确测量高质量角膜图像的角膜厚度,而且针对存在较多噪声和伪影的角膜OCT图像,也可保证角膜厚度测量的准确性。所提角膜厚度测量方法的整体流程如图1所示。

图 1. 所提角膜厚度测量方法流程图

Fig. 1. Flow chart of proposed measurement method of corneal thickness

下载图片 查看所有图片

2.1 图像预处理

图2显示了受噪声和不同伪影源影响的角膜图像。为了消除它们对角膜厚度自动测量的干扰,同时避免边缘模糊,先后采用自适应中值滤波和伪影去除方法对图像进行预处理。

图 2. 原始角膜B扫描图像中存在的噪声和伪影

Fig. 2. Noise and artifacts in original cornea B scan image

下载图片 查看所有图片

2.1.1 自适应中值滤波

在噪声密度较小的情况下(噪声的出现概率小于0.2),使用中值滤波可很好得去除噪声。但是当噪声出现的概率比较大时,中值滤波的效果将受到明显影响。针对存在较多噪声和伪影的角膜OCT图像,采用自适应中值滤波方法来减小图像中噪声的影响。

自适应中值滤波分为以下两个过程:1)判断当前窗口的中值点是否是噪声;2)判断当前窗口的中心位置像素是否是噪声。

如果步骤1)判断中值点不是噪声点,这时跳转到步骤2)进行测试,判断当前窗口的中心位置像素是否为噪声。如果不是噪声,滤波器输出中心位置像素的像素值;如果是噪声,此时输出中值点的像素值。

如果步骤1)判断中值点是噪声点,这时需要增大滤波器的窗口尺寸,在一个更大的范围内寻找一个非噪声的中值,直到找到一个非噪声中值,跳转到步骤2),或者当窗口的尺寸达到了最大值,返回找到的中值。

自适应中值滤波可以处理具有更大概率的脉冲噪声,并且在平滑非脉冲噪声时尽量保留细节,这是传统中值滤波器难以实现的。

2.1.2 去除水平伪影

图2所示的水平伪影是中心饱和伪影与频域OCT直流减除算法作用的结果,被水平伪影影响的行比其他行具有更高的灰度平均强度。目前常用的方法是通过降低水平方向直流信号来消除水平伪影^[1]。该方法先计算出每一行的灰度平均强度,然后再从图像行中的每个像素减去该平均强度。这种方法虽然有效减小了水平伪影,但也降低了图像的对比度和信噪比,不利于角膜图像的边缘提取。

霍夫变换是图像处理中的一种特征提取技术,它利用一种投票算法来检测具有特定形状的物体。该过程在一个参数空间中通过计算累加结果的局部最大值,得到一个符合特定形状的集合作为霍夫变换的结果。在此,采用一种基于霍夫变换的水平伪影去除方法,该处理过程主要分为以下两个步骤:1)利用霍夫变换检测角膜OCT图像中的水平伪影;2) 将图像中水平伪影经过的像素点灰度值置0。

相比于通过降低水平方向直流信号来消除水平伪影的方法,该方法可以在不改变图像对比度和信噪比的情况下,有效地去除图像中的水平伪影。

2.1.3 中心饱和伪影的抑制

中心饱和伪影的关键特征是选定数量的相邻垂直列(A扫描)的强度相对较高。通过搜索相邻垂直列的平均强度的突变,可以确定中心饱和伪影的左右边界。图3显示了图2中A扫描的平均强度曲线,图中两条竖直的黑色虚线表示中央伪影的左右边界,绿色的竖线表示三个相等宽度区域的分界线。先将角膜图像分解成三个相等宽度的区域(图3中的I、II、III),其中只有中心区域包含中心饱和伪影。

图 3. A-扫描平均强度曲线

Fig. 3. A-scan averaged intensity curve

下载图片 查看所有图片

为了确定区域Ⅱ内中心饱和伪影的位置,首先通过计算区域I和区域III垂直列的平均强度来估计没有中心饱和伪影垂直列的平均强度。假设角膜图像的大小为m×n,m为每个A扫描的采样点数,n为A扫描数。每个像素点的灰度值为I(i,j),i为像素点的行坐标,j为像素点的列坐标。区域I所有垂直列的列平均强度为μ₁,区域III所有垂直列的列平均强度为μ₂,则没有中心饱和伪影垂直列的列平均强度μ分别为

μ1=∑i=1m∑j=1n/3I(i,j)mn/3,(1)μ2=∑i=1m∑j=n/3+1nI(i,j)mn/3,(2)μ=μ1+μ22。(3)

在区域Ⅱ中搜索平均强度大于3μ/2的最左垂直列L和最右垂直列R,作为中心饱和伪影的左右边界。在中心饱和伪影左右边界区域内,中心饱和伪影与角膜图像重叠部分的平均强度相对较高。在[L,R]垂直列范围内,计算图像相邻行的平均强度和整个区域的平均强度。相邻行的平均强度曲线如图4所示。

图 4. 相邻行平均强度曲线

Fig. 4. Average intensity curves of adjacent rows

下载图片 查看所有图片

假设第i行的平均强度为P(i),将整个区域的平均强度设定为相邻行平均强度阈值λ_th,搜索平均强度大于λ_th的水平行。根据相邻行的平均强度P(i)是否大于λ_th,每行设置不同的权重,可表示为

ω(i)=0,P(i)<λth1,P(i)≥λth。(4)

然后在[L,R]范围内,将每行的权重ω(i)与对应行每个像素灰度值I(i,j)相乘得到新的像素灰度值I₁(i,j):

I1(i,j)=ω(i)×I(i,j)。(5)

利用该方法对角膜OCT图像进行预处理后的效果如图5所示。图5表明该方法极大地抑制了图像的中心饱和伪影,有利于后续角膜上下边缘轮廓点的提取。

图 5. 角膜B扫描OCT图像经过预处理后的效果。(a)原始角膜B扫描图像;(b)自适应中值滤波后的角膜B扫描图像;(c)去除水平伪影后的角膜B扫描图像;(d)抑制中央伪影的角膜B扫描图像

Fig. 5. Effect of corneal B-scan OCT images after pre-treatment. (a) Original corneal B-scan image; (b) corneal B-scan image with adaptive median-filtering; (c) corneal B-scan image with horizontal artifact removal; (d) corneal B-scan image with suppressed central artifact

下载图片 查看所有图片

2.2 边缘轮廓点提取

采用边界跟踪算法可来提取角膜上下边缘轮廓点,能够准确地提取出角膜的上下边缘轮廓点。在利用边缘跟踪算法获取角膜轮廓点之前,需要将预处理之后的角膜图像进行二值化处理。

2.2.1 图像二值化

图像的二值化有利于图像的进一步处理,使图像变得简单,数据量减少,凸显出感兴趣的目标轮廓。结合坎尼(Canny)边缘检测算子进行图像二值化,能够对低对比度图像和目标像素灰度不均匀的图像较好地进行二值化。具体的处理过程分两个步骤:1)利用Canny算子检测图像边缘,得到角膜边缘图像;2)对角膜边缘图像进行闭运算(先膨胀后腐蚀)处理,得到二值化图像。

2.2.2 边界跟踪算法

边界跟踪算法可提取边界轮廓一系列坐标点或链码,该边界是指一个1像素连通域和0像素连通域之间的边界。二值图像最上面一行、最下面一行、最左边一列和最右边一列组成了图像的边框。灰度值为0和1的像素分别被称为0像素和1像素,假设0像素填充了二值图像的边框,在处理过程中可以给像素点赋任何数值。位于第i行和第j列的像素点用(i,j)表示,采用以下坐标系:行序号i从上到下递增,列序号j从左到右递增。图像像素点(i,j)处的灰度值f_i,j定义为

F={fi,j}。(6)

逐行扫描一幅输入的二值图像,找到一个点(i,j)满足边界跟踪起始点的条件则终止扫描,该点是一个外边界[图6(a)]或是一个孔边界[图6(b)]。如果这个点(i,j)同时满足以上两个条件,(i,j)必须被当作是外边界的起始点。

图 6. 边界条件。(a)外边界;(b)孔边界

Fig. 6. Boundary conditions. (a) Outer boundary condition; (b) hole boundary condition

下载图片 查看所有图片

在此算法中,邻域指的是4-邻域。假设输入图像为F,设置最初的D_NB=1(F的边框形成了一个特殊的孔边界,分配序列数1;D_NB代表当前边界的序列数)。每当开始扫描新的一行,重置D_LNB(上一个遇到的边界序列数)。(i₃,j₃)代表当前点,(i₂,j₂)代表找到的当前点的上一个点,(i₄,j₄)表示在点的邻域中最新找到的非零点,←表示赋值。逐行扫描图像,对每一个满足f_i,j≠0的像素点执行以下步骤。

1) 选择下面的一条:

I如果f_i,j=1且f_i,j-1=0,则确定像素(i,j)是一个外边界跟踪起始点,D_NB自增,并且(i₂,j₂)←(i,j-1);

II如果f_i,j≥1且f_i,j-1=0,则确定像素(i,j)是一个孔边界起始点,D_NB自增,(i₂,j₂)←(i,j+1),当f_i,j>1时,D_LNB←f_i,j;

III否则转到步骤4)。

2) 根据最新找到的边界类型以及带有序列号D_LNB的边界类型,决定当前边界的母边界。

3) 从起始点(i,j)开始,跟踪已检测到的边界,这一步由以下步骤I到步骤V实现。

I从(i₂,j₂)开始,顺时针查找(i,j)的邻域,寻找一个非零点,第一个找到的非零点为(i₁,j₁)。如果没有找到非零点,令f_i,j=-D_NB,转到步骤4)。

II(i₂,j₂)←(i₁,j₁)且(i₃,j₃)←(i,j)。

III寻找范围在当前点(i₃,j₃)的邻域中,方向为逆时针方向,出发点为(i₂,j₂)的下一个元素,将查找的一个非零点设为(i₄,j₄)。

IV根据下面的规则修改点(i₃,j₃)的值f_i3,j3:

a)如果步骤III中检查到点(i₃,j₃+1)是一个零像素点,则f_i3,j3←-D_NB;

b)如果步骤III中检查到点(i₃,j₃)不是一个零像素,且f_i3,j3=1,则f_i3,j3←D_NB;

c)否则,不改变f_i3,j3。

V如果(i₄,j₄)=(i,j)且(i₃,j₃)=(i₁,j₁),转到步骤IV。否则,(i₂,j₂)←(i₄,j₄),(i₃,j₃)←(i₄,j₄),且回到步骤III。

4) 如果f_i,j≠1,则D_LNB← fi,j,且从点(i,j+1)开始重新扫描。当扫描到图片右下角,算法终止。

2.3 角膜厚度的自动测量

图像经过预处理过程后,有效减小了角膜图像的噪声和伪影的影响,但图像中依然存在部分噪声和未消除的伪影,影响角膜轮廓点的提取。为了降低残留噪声和伪影的影响,在利用边界跟踪算法获取图像轮廓之后,寻找角膜轮廓曲线的最小外切矩形,作为角膜在整幅干涉图像中的位置边界。在边界范围内获取角膜的轮廓点,用最小二乘法拟合角膜的上下边界。

假设角膜上边缘拟合曲线方程为y₁=a₁x12+b₁x₁+c₁,下边缘拟合曲线方程为y₂=a₂x22+b₂x₂+c₂。对于角膜上边缘任意一点P₁(x₁,y₁),角膜下边缘上与该点横坐标相同的对应点为P₂(x₂,y₂)。频域OCT(FDOCT)的图像通过干涉光谱的实部傅里叶变换而来,形成的图像不但包括样品本身的层析图像,还包括一个关于零相位延迟对称的镜像。在系统成像深度范围内,可利用的采样点数为每个A扫描采样点数的一半。角膜厚度T的计算公式为

T=2(y2-y1)dN,(7)

式中:d为系统成像深度;N为每个A扫描采样点数。

3 实验

3.1 实验条件

在此采用扫频光学相干层析成像系统获取人眼角膜OCT数据。系统光源中心波长为1060 nm,成像深度为3.7 mm,系统横向分辨率为57.0 μm,纵向分辨率为7.5 μm,信噪比为51 dB,灵敏度为101 dB。每个A扫描采样点数为1280,每幅OCT角膜图像B扫描图像上有360个A扫描。由于频域OCT存在复共轭镜像问题,选取OCT图像的一半作为实验数据,即获得的OCT角膜B扫描图像大小为360 pixel×1280 pixel。采集50幅高质量的角膜OCT图像和50幅存在噪声与伪影的角膜OCT图像,在3.4 GHz的CPU、随机存取存储器为4.0 GB的电脑上利用所提的角膜厚度自动测量方法获得人眼角膜厚度。

3.2 实验结果

人眼角膜OCT B扫描图像的曲线拟合结果如图7所示。由图7可知,利用所提方法不仅可对高质量的角膜OCT图像进行曲线拟合,而且可对存在噪声和伪影的角膜OCT图像进行拟合,为准确测量角膜厚度提供了依据。

图 7. 人眼角膜B扫描OCT图像上下边缘的拟合结果。(a)高质量的原始角膜B扫描图像;(b)高质量的角膜B扫描图像的边缘拟合结果;(c)存在噪声和伪影的原始角膜B扫描图像;(d)存在噪声和伪影的角膜B扫描图像的边缘拟合结果

Fig. 7. Upper and lower edge fitting results of human corneal OCT B-scan image. (a) High quality original corneal B-scan image; (b) high quality edge fitting corneal B-scan image; (c) original corneal B-scan image with noise and artifacts; (d) edge fitting corneal B-scan image with noise and artifacts

下载图片 查看所有图片

使用所提方法分别对某健康人50幅高质量角膜OCT图像和50幅存在噪声和伪影的角膜OCT图像进行曲线拟合后,计算角膜沿Y轴方向的平均厚度(T_a)、中央厚度(T_c)和对应的标准方差(σ),结果如表1所示。

从表1中可以看出,两种类型的角膜OCT图像中沿Y轴方向的平均厚度相差0.5 μm,中央厚度相差0.7 μm。噪声和伪影对角膜边缘曲线的拟合以及角膜中央厚度的测量有一定的影响,但沿Y轴方向平均厚度的偏差和中央厚度的偏差都非常小,可以认为测量偏差处于可接受范围内。由此证明,该方法适用于存在噪声和伪影的角膜OCT图像角膜厚度的自动测量。

4 分析与讨论

4.1 不同预处理去噪算法对比

目前常用的图像去噪方法主要有均值滤波、中值滤波和高斯滤波,但这几种方法均存在一定的缺陷。均值滤波无法完全消除图像中的噪声,减弱噪声的效果受限;中值滤波容易导致图像的不连续性;高斯滤波对于抑制服从正态分布的噪声非常有效,但不能很好地去除图像中的散斑噪声。为了解决这几种常用的去噪方法存在的问题,在此采用基于自适应中值滤波的噪声消除算法去除图像中的噪声。图8显示的是未进行去噪预处理、使用均值滤波进行预处理、使用中值滤波进行预处理以及使用自适应中值滤波进行预处理等几种情况下,利用边界跟踪与曲线拟合相结合的方法得到的OCT角膜B扫描图像上下边缘的曲线拟合效果。计算经过不同预处理去噪算法处理后的角膜厚度,并将其与人工手动测量的角膜厚度进行对比,对比结果如表2所示。

图 8. 不同预处理去噪算法得到的人眼角膜B扫描图像上下边缘的拟合效果。(a)未进行预处理去噪;(b)均值滤波预处理;(c)中值滤波预处理;(d)自适应中值滤波预处理

Fig. 8. Fitting effects of upper and lower edges of human corneal B-scan images obtained by different pre-processing denoising algorithms. (a) No preprocessing denoising. (b) mean filter preprocessing; (c) median filter preprocessing; (d) adaptive median filter preprocessing

下载图片 查看所有图片

对图8(a)~(c)中箭头所指区域进行拟合,得到的曲线偏离了角膜B扫描图像的上下边缘,对角膜厚度的准确测量产生了较大的影响。图8(d)中的拟合曲线未偏离角膜B扫描图像的上下边缘,得到的拟合曲线更加平滑、准确,保证了角膜测量的精确性。表2中以人工手动测量的角膜厚度作为基准,与其他预处理去噪算法相比,经过自适应中值滤波处理后的角膜厚度与人工手动测量的角膜厚度偏差最小,并且对应的平均方差也最小。因此,自适应中值滤波的预处理效果更好。

4.2 与基于边缘检测、随机抽样一致性方法的对比

为了评价所提出的角膜厚度自动测量方法,使用基于边缘检测和随机抽样一致性的方法与所提方法自动测量相同的OCT角膜B扫描图像中的角膜厚度,并对最终的测量结果进行对比。为了确保随机抽样一致性算法的正确性,按照2012年舒鹏等^[7]提出的方法完成角膜厚度的自动测量。所提方法与基于边缘检测和随机抽样一致性的方法对相同OCT角膜B扫描图像上下边缘的曲线拟合效果对比如图9所示。

图 9. OCT角膜图像上下边缘的拟合结果。(a)边缘检测和随机抽样一致性方法;(b)所提方法

Fig. 9. Results of upper and lower edge fitting of OCT corneal images. (a) Edge detection and random sampling consistency method; (b) proposed method

下载图片 查看所有图片

由图9可知,所提方法可以准确拟合出角膜上下边缘,而基于边缘检测和随机抽样一致性方法的角膜上下边缘拟合结果存在一定偏差,在图中箭头所指区域拟合曲线偏离角膜上边缘较大。与基于边缘检测和随机抽样一致性方法相比,所提方法拟合角膜上下边缘的精度更高,在角膜图像低信噪比和存在中央伪影的区域也有较好的拟合效果,对噪声的抗干扰能力更强。

分别将两种方法的角膜厚度计算结果与医生手动选取采样点得到的角膜厚度结果进行对比^[25-26],计算得到沿Y轴方向的平均厚度偏差f_mean(A,M)、角膜中央厚度偏差以及对应的标准偏差f_std(A,M),计算结果如表3所示,计算公式为

fmean(A,M)=1uv∑k=1u∑l=1vA(k,l)-M(k,l),(8)fstd(A,M)=1u∑k=1u1v∑l=1vA(k,l)-M(k,l)-fmean(A,M)2,(9)

式中:u为角膜B扫描图像A扫描数;v为角膜B扫描幅数;A_(k,l)和M_(k,l)分别表示第l幅角膜B扫描图像中第k个A扫描上对应的角膜厚度的自动测量结果和人工计算结果。将u置为1,k置为角膜中央所在A-line,计算结果即为角膜中央厚度偏差。

由表3可知,利用所提方法得到的角膜沿Y轴方向的平均厚度偏差、中央厚度偏差和对应的标准方差均小于基于边缘检测和随机抽样一致性方法的相应结果。图9(a)中的图像在低信噪比和存在中央伪影的区域,角膜上边缘的曲线拟合结果存在较大偏离,表3中的数据也验证了图9的分析结果。总体来说,所提方法适应于低信噪比和含高噪声的图像,对存在水平和中央伪影的图像也具有很好的测量结果。

5 结论

提出了一种角膜OCT图像中角膜厚度自动测量方法,该方法采用自适应中值滤波、霍夫直线检测和中央伪影去除的方法对角膜OCT图像进行预处理,有效减小图像中散斑噪声和伪影的影响;然后利用边界跟踪算法和轮廓定位的方法获取角膜上下边缘轮廓点,根据获取的轮廓点并利用最小二乘法拟合出角膜上下边缘曲线,得到拟合曲线方程,最终计算出角膜厚度。与已有的角膜厚度测量方法相比,所提方法可有效减小噪声和伪影对角膜厚度测量的影响,角膜上下边缘曲线拟合精度更高,角膜厚度测量更准确。