全文HTML
--> --> --> -->2.1.角膜高、低信噪比区域划分
为实现对角膜图像的分区处理, 首先需要将其划分为高、低信噪比区域. 假设角膜图像大小为M × N, M为A-scan的采样点数, N为A扫描数, 每个像素点的灰度值为I(i,j), 其中i为像素点的行坐标, j为像素点的列坐标. 利用(1)式计算各列像素间的标准差, 为便于观察, 对(1)式的结果进行平滑滤波, 处理结果如图2(a)所示, 其中横轴信息为像素列坐标, 纵轴信息为各列像素的标准差平滑结果. 计算图像各列像素标准差的均值meanall, 并以meanall为阈值对图像进行高、低信噪比区域的分割. 在图2(a)中搜索纵坐标大于meanall的区域, 最终定位高信噪比区域为Region II, 低信噪比区域为Region I, Region III.图 2 角膜图像高、低信噪比区域划分 (a) 列间标准差平滑结果; (b) 角膜图像高、低信噪比区域划分结果
Figure2. Division between high and low SNR (signal-to-noise ratio) regions of corneal image: (a) Smoothing result of standard deviation between columns; (b) division results of high and low SNR regions of corneal image.
2
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)所示, 其中红色点表示上表面轮廓点, 绿色点表示下表面轮廓点. 初步提取结果存在一定误差, 其主要来源于两部分: 上下表面灰度分布不均造成轮廓点定位结果与角膜实际结构不符; 孤立噪点引起的定位错误.
图 3 高信噪比区域的轮廓提取结果 (a) 轮廓点初步提取结果; (b) 轮廓点精确提取结果
Figure3. Contour extraction results of high SNR region: (a) Preliminary extraction result of contour points; (b) accurate extraction result of contour points.
针对第一项误差, 可依据轮廓点定位结果的位置信息即坐标值大小将其划分为上、下表面; 针对第二项误差, 可根据角膜结构具有局部稳定性[23]的特征, 计算轮廓点中的所有相邻像素的坐标差值, 取其平均值为diff_mean, 根据多次实验结果, 当某一轮廓点与其相邻像素的坐标差值超过3diff_mean即可认为该点为误差点.
经由上述算法的处理结果如图3(b)所示, 可以看出通过上述算法, 高信噪比区域的角膜轮廓得到精确提取, 伪影部分及缺损结构经由其相邻像素的补充也得到了相对准确的定位.
2
2.3.图像增强
单帧角膜图像低信噪比区域与背景信号基本融为一体, 信噪比极低, 无法直接利用灰度分布来获取角膜轮廓. 由于SS-OCT系统成像速度极快, 连续帧角膜图像除少量刚性位移外, 基本可以认为角膜结构不变; 而角膜图像中存在包含散斑噪声[24]在内的多种随机噪声, 因此可以利用对噪声容忍度较高的相位相关[25]算法对连续帧图像进行配准叠加, 帧间叠加过程中利用噪声的随机性对噪声进行抑制同时实现对低信噪比区域角膜轮廓的增强.经实验验证, 对连续5帧角膜图像进行配准叠加即可在保持角膜结构不变的前提下实现角膜图像轮廓的增强. 由于角膜图像直接叠加会使得伪影区域大面积扩增并出现过饱和现象, 造成伪影与角膜轮廓混叠以至于无法区分. 因此首先利用2.1节分别对单帧图像高、低信噪比区域进行划分, 接着利用2.2节算法分别对单帧图像进行高信噪比区域的轮廓定位, 最后对定位结果进行配准叠加, 得到待处理的角膜增强图像, 实验结果表明, 经配准叠加后的角膜图像低信噪比区域的信噪比从原来的单帧图像的信噪比更高.
2
2.4.图像背景噪声的滤除
由于角膜增强图像中大部分像素为背景信息, 因此可通过OTSU算法提取出角膜轮廓信息, 提取结果如图4(a)所示. 由于上述滤波结果中仍存在孤立噪声点, 考虑到需要保持图像的纵向轮廓信息, 而角膜内部的细节信息在轮廓提取时并不重要. 因此对图4(a)结果进行横向中值滤波, 以实现对角膜图像噪声的抑制, 滤波结果如图4(b)所示.图 4 角膜整体轮廓提取过程 (a) OTSU算法处理结果; (b) 中值滤波处理结果
Figure4. Extraction process of the overall cornea contour: (a) Processing result by OTSU algorithm; (b) processing result by median filtering.
2
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, 则针对角膜上、下表面, 分别进行以下判断:
图 5 低信噪比区域角膜轮廓建模及提取结果 (a)低信噪比区域表面轮廓点建模结果; (b)低信噪比区域轮廓提取结果
Figure5. Modeling and extraction results of corneal contour in low SNR region: (a) Modeling results of surface contour points in low SNR region; (b) contour extraction result of low SNR region.
(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
2.6.角膜完整轮廓提取
经过2.1—2.4节提取的高、低信噪比区域的轮廓点(图6(a))并非平滑连续的角膜结构, 因此需要对提取的轮廓点进行全局多项式拟合, 通过整合高、低信噪比轮廓点的位置信息, 以得到完整的角膜轮廓. 拟合结果如图6(b)所示.图 6 角膜完整轮廓的提取过程 (a)角膜上下表面轮廓点提取结果; (b) 角膜轮廓拟合结果
Figure6. Extraction process of the complete cornea contour: (a) Extraction results of the contour points in the upper and lower cornea surfaces; (b) fitting result of the cornea contour
在现有算法中, 只有文献[14]能够同时实现对含伪影及低信噪比区域的角膜图像轮廓的提取, 但其对低信噪比区域的处理采用传统的二项式拟合[14], 所获取的弱边缘轮廓精度较低. 而本文算法利用了低信噪比区域的真实信息引导轮廓结构的提取, 可增强提取精度与可靠性.
为验证本文算法能否实现对低质量角膜图像进行精确的轮廓提取, 分别利用文献[14]中的算法与本文算法对获取的10组角膜图像进行处理, 两种算法在实验平台上的平均数据处理时间分别为9.6 s和8.9 s.
处理的其中一组结果如图7(a)所示, 其他组实验与该结果类似. 图中的红色和玫红色曲线为本文算法的提取结果, 蓝色和绿色曲线为文献[14]提取结果, 局部图A为角膜上表面提取结果, 局部图B为角膜下表面提取结果, 局部图C为含伪影部分的高信噪比区域提取结果; 分别计算两种算法获取的角膜各位置的平均厚度, 计算结果如图7(b)所示, 图中蓝色曲线表示本文算法获取的角膜厚度值与Pstd间的偏差随位置信息的变化情况, 绿色曲线为文献[14]获取的角膜厚度值与Pstd间的偏差随位置信息的变化情况, 其中, 红色直线为角膜实际厚度Pstd; 分别计算两种算法在高、低信噪比区域提取结果与Pstd间的平均误差值与偏差率, 结果如表1所示.
图 7 两种算法效果对比 (a)两种算法轮廓提取结果; (b) 两种算法角膜厚度平均计算结果
Figure7. Comparison of the effects of the two algorithms: (a) Results of the contour extraction of the two algorithms; (b) results of the corneal thickness calculated by the two algorithms.
评价指标 | 文献[14] | 本文算法 | 精度提高值/% |
平均误差值/pixels | 2.8 | 0.9 | 4.9 |
偏差率/% | 7.2 | 2.3 |
表1两种算法对角膜轮廓平均提取精度对比
Table1.Comparison of the accuracy of two algorithmsfor contour extraction of high and low SNR regions.
由图7(a)的小图A, B及图7(b)可以看出, 在低信噪比区域, 本文算法提取的角膜轮廓更加贴近角膜上下表面. 由图7(a)的C可以看出, 在高信噪比区域, 本文算法与文献[14]的提取结果相近, 两种算法相较于角膜的实际标准尺寸的平均测量偏差分别为7.2%和2.3%. 本文算法获取的角膜厚度更加接近角膜实际厚度, 平均提取精度较文献[14]的提高了4.9%.