1.   引　言

扫频光学相干层析成像(SS-OCT)是一种基于低相干干涉测量的医学成像技术^[1-3], 具有非侵入、分辨率高、成像速度快的优势, 近年来, 已成为临床眼科检查角膜的一种标准手段^[4,5]. SS-OCT系统获取的角膜轮廓包含角膜曲率及角膜厚度等信息^[6], 这些信息是角膜异常定量分析及眼科临床手术的重要依据^[7-9]. 由于角膜图像中常含有低信噪比区域与中心饱和伪影, 同时存在轮廓信息缺失及模糊现象^[10], 而完整精准的角膜轮廓能够提供更加可靠的病理信息, 因此, 实现对低质量的角膜图像轮廓的精确提取十分必要.

传统的角膜轮廓提取方法采用手动追踪, 但该过程冗长耗时且易出现人为误差. 为克服上述困难, 多种自动提取算法被提出. Li等^[11,12]提出结合快速活动轮廓和二阶多项式拟合的方法, 但该方法未考虑中央伪影且易受噪声干扰; Shu和Sun^[13]提出基于边缘检测和随机抽样一致的方法, 虽然排除了中央伪影的影响, 但未考虑角膜弱边缘的存在; LaRocca等^[14]提出鲁棒性更强的基于图论和动态规划^[15]的方法, 虽然能够实现对同时包含伪影及低信噪比区域的角膜轮廓的提取, 但其对低信噪比区域的定位采用Li等^[11,12]提出的二阶多项式近似的假设, 定位结果并不精确; Santos等^[16]通过深度学习的方法实现对角膜的分割, 但是该方法需要大量临床数据及专业处理平台.

本文提出一种针对SS-OCT系统获取的低质量角膜图像的轮廓自动提取算法, 先根据图像标准差将角膜图像划分为高、低信噪比区域并分别进行轮廓提取, 再通过全局多项式拟合实现区域间的轮廓合并与重整.

2.   算法原理

SS-OCT成像系统获取的角膜图像因其远心扫描模式造成图像中央区域信噪比较高, 轮廓线条明显, 但含有伪影区域; 角膜两侧信噪比较低, 角膜轮廓较为模糊, 无法单纯利用图像梯度信息进行提取, 同时SS-OCT系统的振镜组高速扫描模式也造成角膜图像轮廓及对比度并非连续均匀的现象. 因此, 针对包含伪影及部分结构缺失的高信噪比区域, 可根据像素的灰度分布来定位峰值点作为角膜轮廓; 针对轮廓信息基本缺失的低信噪比区域, 可通过帧间叠加对角膜图像进行弱边缘增强从而提供部分低信噪比区域的轮廓信息, 同时利用该信息引导角膜轮廓局部直线拟合结果, 使其更加贴近角膜真实轮廓; 针对角膜轮廓非连续均匀的问题, 可通过全局多项式拟合实现对高、低信噪比区域角膜轮廓的整合及提取. 角膜图像轮廓算法流程如图1所示.

首先根据图像标准差分布将图像划分为高、低信噪比区域; 之后对角膜图像高信噪比区域进行轮廓提取; 接着利用连续多帧图像叠加进行图像增强, 并通过最大类间方差法(OTSU)^[17]及中值滤波算法实现对图像背景噪声的滤除; 之后通过权衡局部直线拟合结果及参考轮廓点实现对角膜低信噪比区域的轮廓提取; 最后, 对高、低信噪比区域的所有轮廓点进行全局多项式拟合, 实现对角膜轮廓的自动提取.

2.1   角膜高、低信噪比区域划分

为实现对角膜图像的分区处理, 首先需要将其划分为高、低信噪比区域. 假设角膜图像大小为M × N, M为A-scan的采样点数, N为A扫描数, 每个像素点的灰度值为I(i,j), 其中i为像素点的行坐标, j为像素点的列坐标. 利用(1)式计算各列像素间的标准差, 为便于观察, 对(1)式的结果进行平滑滤波, 处理结果如图2(a)所示, 其中横轴信息为像素列坐标, 纵轴信息为各列像素的标准差平滑结果. 计算图像各列像素标准差的均值mean_all, 并以mean_all为阈值对图像进行高、低信噪比区域的分割. 在图2(a)中搜索纵坐标大于mean_all的区域, 最终定位高信噪比区域为Region II, 低信噪比区域为Region I, Region III.

$$\begin{split}{\rm{std}}(j) = \; & \sqrt {\frac1{M}{{\displaystyle\sum\limits_{i = 1}^M {{{(I(i,j) - {{\bar I}_j})}^2}} }}},\\ & j = 1,2 \cdots N,~~i = 1,2 \cdots M.\end{split}$$

(1)

2.2   高信噪比区域的轮廓提取

由于高信噪比区域轮廓线条明显, 其灰度明显高于背景信号, 一般表现为单列像素中的峰值像素点. 因此对于高信噪比区域的轮廓, 可根据其单列像素的灰度分布进行峰值点提取从而获取角膜轮廓. 由于角膜轮廓部分结构缺失及伪影区域的存在, 导致单列像素峰值点提取过程中存在以下几种情况:

1)角膜上、下表面均不含伪影且不存在结构缺失;

2)角膜仅有一表面含有伪影;

3)角膜上、下表面均含有伪影;

4)角膜某一表面或上、下表面边缘结构缺失.

为解决不同情况下对轮廓点的定位, 首先设定图像灰度均值threshold为判断像素点是否为轮廓点的阈值, 并依据人眼常见病变如高眼压症^[18]、远视眼^[19]、近视眼^[20,21]、青光眼^[22]等以及人眼正常生理性变化^[23], 设定角膜上下表面轮廓点间距范围为CT_range, 伪影区域范围为Artifact_range.

针对情况1), 可直接搜索单列像素中大于threshold的峰值点及与峰值点间距大于CT_range的次峰值点作为角膜轮廓点. 针对情况2), 在搜索大于threshold的峰值点时, 若同时存在多个像素点灰度相同或极为接近, 且集中在Artifact_range范围内, 说明该像素点集处于伪影区域, 则选取像素点集的中值点(针对灰度相同)或峰值点(针对灰度接近)作为轮廓点, 其余点认定为伪影点, 后续不予处理, 并继续在该列残余像素中搜索与已定位的轮廓点间距大于CT_range的次峰值点作为另一表面轮廓点. 针对情况3), 在搜索大于threshold的峰值点时, 若同时存在多个像素点灰度相同, 且像素点中存在相邻点对的间距超过CT_range时, 则说明上、下表面轮廓点均存在伪影且伪影像素点与峰值点混叠, 因此以相邻点对为界将像素点集分为两部分, 并选取各部分中值点作为轮廓点, 以保证伪影区域的轮廓点尽可能接近真实轮廓. 针对情况4), 当角膜结构缺失时, 说明该列像素中不存在大于threshold的峰值点或存在峰值点但不存在与其间距超过CT_range的次峰值点, 针对这种情况, 可依据后续的局部多项式拟合算法进行缺失结构的填充.

根据上述方法对角膜进行初步轮廓提取的结果如图3(a)所示, 其中红色点表示上表面轮廓点, 绿色点表示下表面轮廓点. 初步提取结果存在一定误差, 其主要来源于两部分: 上下表面灰度分布不均造成轮廓点定位结果与角膜实际结构不符; 孤立噪点引起的定位错误.

针对第一项误差, 可依据轮廓点定位结果的位置信息即坐标值大小将其划分为上、下表面; 针对第二项误差, 可根据角膜结构具有局部稳定性^[23]的特征, 计算轮廓点中的所有相邻像素的坐标差值, 取其平均值为diff_mean, 根据多次实验结果, 当某一轮廓点与其相邻像素的坐标差值超过3diff_mean即可认为该点为误差点.

经由上述算法的处理结果如图3(b)所示, 可以看出通过上述算法, 高信噪比区域的角膜轮廓得到精确提取, 伪影部分及缺损结构经由其相邻像素的补充也得到了相对准确的定位.

2.3   图像增强

单帧角膜图像低信噪比区域与背景信号基本融为一体, 信噪比极低, 无法直接利用灰度分布来获取角膜轮廓. 由于SS-OCT系统成像速度极快, 连续帧角膜图像除少量刚性位移外, 基本可以认为角膜结构不变; 而角膜图像中存在包含散斑噪声^[24]在内的多种随机噪声, 因此可以利用对噪声容忍度较高的相位相关^[25]算法对连续帧图像进行配准叠加, 帧间叠加过程中利用噪声的随机性对噪声进行抑制同时实现对低信噪比区域角膜轮廓的增强.

经实验验证, 对连续5帧角膜图像进行配准叠加即可在保持角膜结构不变的前提下实现角膜图像轮廓的增强. 由于角膜图像直接叠加会使得伪影区域大面积扩增并出现过饱和现象, 造成伪影与角膜轮廓混叠以至于无法区分. 因此首先利用2.1节分别对单帧图像高、低信噪比区域进行划分, 接着利用2.2节算法分别对单帧图像进行高信噪比区域的轮廓定位, 最后对定位结果进行配准叠加, 得到待处理的角膜增强图像, 实验结果表明, 经配准叠加后的角膜图像低信噪比区域的信噪比从原来的单帧图像的信噪比更高.

2.4   图像背景噪声的滤除

由于角膜增强图像中大部分像素为背景信息, 因此可通过OTSU算法提取出角膜轮廓信息, 提取结果如图4(a)所示. 由于上述滤波结果中仍存在孤立噪声点, 考虑到需要保持图像的纵向轮廓信息, 而角膜内部的细节信息在轮廓提取时并不重要. 因此对图4(a)结果进行横向中值滤波, 以实现对角膜图像噪声的抑制, 滤波结果如图4(b)所示.

2.5   低信噪比区域的轮廓提取

以高、低信噪比区域的边界为起始搜索位置进行轮廓点的搜索, 以待定位点的前向多个已定位轮廓点为基础进行局部直线拟合, 并以此预测该待定位点的坐标, 同时权衡2.3节及2.4节提供的部分低信噪比区域参考轮廓结构信息与该预测值的优劣, 确定最佳匹配轮廓点. 设定最佳匹配轮廓点坐标为res, 各像素点前向搜索范围为range, 以预测值为原点的角膜轮廓上下搜索范围为2tolerance, 经实验验证, 本文选取的range = CT_range/4, tolerance = range/5. 具体实现方法如下.

1)对第j列前向range列内的已定位轮廓点进行直线拟合, 得到预测直线polyfit, 由polyfit得到第j列的轮廓点预测值pred(j).

2)对低信噪比区域表面轮廓点进行建模(如图5(a)所示), 在第j列角膜各表面tolerance范围内搜索非零像素点即参考轮廓点nonzero, 若存在nonzero, 则针对角膜上、下表面, 分别进行以下判断:

(a)对于上表面, 若nonzero点坐标大于等于pred(j), 即针对参考轮廓点, 预测值更接近上表面, 则pred(j)为最佳匹配轮廓点(见第11, 10列); 若nonzero点坐标小于pred(j), 则nonzero为最佳匹配轮廓点(见第8列);

(b)对于下表面, 若nonzero点坐标小于等于pred(j), 即针对参考轮廓点, 预测值更接近下表面, 则pred(j)为最佳匹配轮廓点(见第8, 10列); 若nonzero点坐标大于pred(j), 则nonzero为最佳匹配轮廓点(见第11列).

3)若搜索范围内不存在nonzero, 则obs(j)为最佳匹配轮廓点(见第9列).

通过上述方法实现对低信噪比区域的轮廓提取, 最终得到的角膜边缘如图5(b)所示.

2.6   角膜完整轮廓提取

经过2.1—2.4节提取的高、低信噪比区域的轮廓点(图6(a))并非平滑连续的角膜结构, 因此需要对提取的轮廓点进行全局多项式拟合, 通过整合高、低信噪比轮廓点的位置信息, 以得到完整的角膜轮廓. 拟合结果如图6(b)所示.

3.   实验及结果

为验证本文算法的正确性, 首先搭建一套SS-OCT成像系统对角膜进行成像. 该系统选用中心波长为1060 nm, 带宽为30 nm, 频率为3 kHz, 相干长度为69 mm的扫频光源(AXSUN-1060), 其在组织内的理论纵向分辨率为15 μm. 样品选用光学眼模型, 其角膜折射率为1.376, 角膜厚度为550 μm, 符合正常成人角膜结构参数. 对该样品连续采集5帧图像作为一组, 共计采集10组图像, 采集的每幅图像均包含820条A-scan, 每条A-scan包含350个像素点, 根据系统分辨率计算出角膜实际厚度对应的像素点数为P_std. 本文采用Intel Core i5-3337U处理器, 4G内存, 以Windows x64系统下的MATLAB 2017b为平台进行实验.

在现有算法中, 只有文献[14]能够同时实现对含伪影及低信噪比区域的角膜图像轮廓的提取, 但其对低信噪比区域的处理采用传统的二项式拟合^[14], 所获取的弱边缘轮廓精度较低. 而本文算法利用了低信噪比区域的真实信息引导轮廓结构的提取, 可增强提取精度与可靠性.

为验证本文算法能否实现对低质量角膜图像进行精确的轮廓提取, 分别利用文献[14]中的算法与本文算法对获取的10组角膜图像进行处理, 两种算法在实验平台上的平均数据处理时间分别为9.6 s和8.9 s.

处理的其中一组结果如图7(a)所示, 其他组实验与该结果类似. 图中的红色和玫红色曲线为本文算法的提取结果, 蓝色和绿色曲线为文献[14]提取结果, 局部图A为角膜上表面提取结果, 局部图B为角膜下表面提取结果, 局部图C为含伪影部分的高信噪比区域提取结果; 分别计算两种算法获取的角膜各位置的平均厚度, 计算结果如图7(b)所示, 图中蓝色曲线表示本文算法获取的角膜厚度值与P_std间的偏差随位置信息的变化情况, 绿色曲线为文献[14]获取的角膜厚度值与P_std间的偏差随位置信息的变化情况, 其中, 红色直线为角膜实际厚度P_std; 分别计算两种算法在高、低信噪比区域提取结果与P_std间的平均误差值与偏差率, 结果如表1所示.

由图7(a)的小图A, B及图7(b)可以看出, 在低信噪比区域, 本文算法提取的角膜轮廓更加贴近角膜上下表面. 由图7(a)的C可以看出, 在高信噪比区域, 本文算法与文献[14]的提取结果相近, 两种算法相较于角膜的实际标准尺寸的平均测量偏差分别为7.2%和2.3%. 本文算法获取的角膜厚度更加接近角膜实际厚度, 平均提取精度较文献[14]的提高了4.9%.

4.   总　结

本文提出了一种针对低质量的SS-OCT角膜图像的轮廓自动提取算法, 无需手动追踪角膜轮廓, 直接通过对高、低信噪比区域分别处理即可提取出低质量角膜图像的角膜轮廓信息. 该算法通过帧间叠加实现了图像增强, 继而引入了角膜的实际轮廓信息作为权衡因子, 增强了算法的鲁棒性. 与目前仅有的能够处理既存在伪影又存在低信噪比区域图像的角膜轮廓提取算法相比, 本文算法对低信噪比区域角膜轮廓的提取精度更高, 能够为后续低质量角膜图像轮廓的精确提取提供重要依据.