留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

扫频光学相干层析视网膜图像配准去噪算法

蔡怀宇 韩晓艳 娄世良 汪毅 陈文光 陈晓冬

蔡怀宇, 韩晓艳, 娄世良, 汪毅, 陈文光, 陈晓冬. 扫频光学相干层析视网膜图像配准去噪算法[J]. 中国光学. doi: 10.37188/CO.2020-0130
引用本文: 蔡怀宇, 韩晓艳, 娄世良, 汪毅, 陈文光, 陈晓冬. 扫频光学相干层析视网膜图像配准去噪算法[J]. 中国光学. doi: 10.37188/CO.2020-0130
CAI Huai-yu, HAN Xiao-yan, LOU Shi-liang, WANG Yi, CHEN Wen-guang, CHEN Xiao-dong. Speckle noise reduction in swept-source optical coherence tomography by retinal image registration[J]. Chinese Optics. doi: 10.37188/CO.2020-0130
Citation: CAI Huai-yu, HAN Xiao-yan, LOU Shi-liang, WANG Yi, CHEN Wen-guang, CHEN Xiao-dong. Speckle noise reduction in swept-source optical coherence tomography by retinal image registration[J]. Chinese Optics. doi: 10.37188/CO.2020-0130

扫频光学相干层析视网膜图像配准去噪算法

doi: 10.37188/CO.2020-0130
基金项目: 国家重点研发计划(No.2017YFC0109901);天津市自然科学基金项目(No.15JCQNJC14200)
详细信息
    作者简介:

    蔡怀宇(1965—),女,湖南涟源人,博士,教授,硕士生导师,1991年、2000年于天津大学分别获得硕士、博士学位,主要从事信息光学、光电技术及仪器和图像处理等方面的研究。E-mail:hycai@tju.edu.cn

    韩晓艳(1996—),女,内蒙古乌兰察布人,天津大学精密仪器与光电工程技术学院硕士研究生,2018年于南京理工大学获得学士学位,主要从事光学相干层析成像方面的研究。E-mail:xyhan@tju.edu.cn

  • 中图分类号: TN247

Speckle noise reduction in swept-source optical coherence tomography by retinal image registration

Funds: Supported by National Key R&D Program of China (No. 2017YFC0109901); Natural Science Foundation Project of Tianjin (No. 15JCQNJC14200)
More Information
  • 摘要: 多帧叠加平均处理是去除扫频光学相干层析系统的散斑噪声、获得较为清晰结构信息的有效方法,但眼睛的震颤、漂移、微眼跳等生理特性和系统光路特性会使图像之间存在几何变换,导致叠加效果不佳、结构稳定性差,为此本文提出一种基于灰度分布信息和目标几何信息相结合的配准算法。该方法根据图像平均灰度分布提取包含目标信息的感兴趣区域,通过相位相关算法和基于分段拟合的灰度投影算法的双重作用校正图像的平移变换;通过拟合视网膜上边界作为特征点迭代确定最佳旋转参数,并再次重新估计平移参数,实现图像的刚性配准;最后通过轴向扫描一对一映射法以能量函数为约束条件实现图像的非刚性配准。对活体兔眼进行实验,结果表明,本文算法配准后的叠加图像边界清晰,结构信息增强,信噪比和对比度平均有效提高一倍多。本算法适用于强噪声视网膜B-Scans图像的配准,能满足多种类型OCT系统的叠加成像需要,具有较高的鲁棒性和图像配准精度。
  • 图  1  OCT序列图像几何变换原因分析(a)眼睛轴向运动(b)眼睛垂轴运动(c)眼睛旋转运动(d)多角度光线入射

    Figure  1.  Analysis of reasons for geometric transformation of OCT sequence images (a) Axial eye movement (b) Vertical eye movement (c) Rotational eye movement (d) Multi-angle light incidence

    图  2  配准算法流程图

    Figure  2.  Flow chart of the registration algorithm

    图  3  预处理过程(a)原始视网膜B-Scan图像(b)列平均灰度分布图(c)ROI图像

    Figure  3.  Pretreatment process (a) Original B-Scan retinal image (b) Column average gray distribution map (c) ROI image

    图  4  分段拟合处理过程(a)曲线分割结果(黑色虚线为分割线)(b)原始曲线段(最佳连接点A、黑色直线为包络线)(c)处理后的曲线段

    Figure  4.  Fitting process of curve segments (a) Results of curve segmentation (the dotted black line is the segmentation line) (b) Original curve segment (optimal connection point A, the straight black line is the envelope) (c) Curve segment after processing

    图  5  刚性配准结果(a)上边界提取过程(b)平移配准后的叠加图(c)旋转配准后的叠加图(绿色为参考图像,紫色为目标图像)

    Figure  5.  Results of rigid registration (a) Upper boundary extraction process (b) Averaging image after translation registration (c) Averaging image after rotation registration (Green is the reference image, purple is the target image)

    图  6  非刚性配准结果(a)刚性配准后的叠加图及局部放大图(b)非刚性配准后的叠加图及局部放大图(绿色为参考图像,紫色为目标图像)

    Figure  6.  Results of non-rigid registration (a) Average image obtained after rigid registration (b) Average image obtained after non-rigid registration (The red dashed box represents an enlarged partial view; Green is the reference image, and purple is the target image)

    图  7  本算法实验结果图(a)原始单帧图(b)配准后的叠加图

    Figure  7.  The experimental results of the algorithm (a) Original single frame image (b) Averaging image after registration

    图  8  SNR和CNR随叠加帧数递增的曲线变化图(a)区域选择(b)SNR(c)CNR

    Figure  8.  SNR and CNR curve change graph with an increasing number of multiple frames (a) Area selection (b) SNR (c) CNR

    表  1  配准算法对比结果

    Table  1.   Comparison of the registration algorithms

    比较对象SNR提高倍数CNR提高倍数
    本文方法/文献[18]方法1.121.19
    本文方法/文献[21]方法1.501.48
    下载: 导出CSV
  • [1] SAHYOUN C C, SUBHASH H M, PERU D, et al. An experimental review of optical coherence tomography systems for noninvasive assessment of hard dental tissues[J]. Caries Research, 2019, 54(1): 43-54.
    [2] PASAOGLU I, SATANA B, ALTAN C, et al. Lamina cribrosa surface position in idiopathic intracranial hypertension with swept-source optical coherence tomography[J]. Indian Journal of Ophthalmology, 2019, 67(7): 1085-1088. doi:  10.4103/ijo.IJO_1736_18
    [3] PARK J H, YOO C, JUNG J H, et al. The association between prelaminar tissue thickness and peripapillary choroidal thickness in untreated normal-tension glaucoma patients[J]. Medicine, 2019, 98(1): e14044. doi:  10.1097/MD.0000000000014044
    [4] 汪毅, 刘珊珊, 张玮茜, 等. 扫频光学相干层析角膜图像轮廓自动提取算法[J]. 物理学报,2019,68(20):204201. doi:  10.7498/aps.68.20190731

    WANG Y, LIU SH SH, ZHANG W Q, et al. Automatic contour extraction algorithm for swept-source optical coherence tomography cornea image[J]. Acta Physica Sinica, 2019, 68(20): 204201. (in Chinese) doi:  10.7498/aps.68.20190731
    [5] VIRA J, MARCHESE A, SINGH R B, et al. Swept-source optical coherence tomography imaging of the retinochoroid and beyond[J]. Expert Review of Medical Devices, 2020, 17(5): 413-426. doi:  10.1080/17434440.2020.1755256
    [6] LOU SH L, CHEN X D, HAN X Y, et al. Fast retinal segmentation based on the wave algorithm[J]. IEEE Access, 2020, 8: 53678-53686. doi:  10.1109/ACCESS.2020.2981206
    [7] SCHMITT J M, XIANG S H, YUNG K M, et al. Speckle in optical coherence tomography[J]. Journal of Biomedical Optics, 1999, 4(1): 95-105. doi:  10.1117/1.429925
    [8] 蔡怀宇, 张玮茜, 陈晓冬, 等. 眼科光学相干层析成像的图像处理方法[J]. 中国光学,2019,12(4):731-740. doi:  10.3788/co.20191204.0731

    CAI H Y, ZHANG W Q, CHEN X D, et al. Image processing method for ophthalmic optical coherence tomography[J]. Chinese Optics, 2019, 12(4): 731-740. (in Chinese) doi:  10.3788/co.20191204.0731
    [9] JORGENSEN T M, THOMADSEN J, CHRISTENSEN U, et al. Enhancing the signal-to-noise ratio in ophthalmic optical coherence tomography by image registration—method and clinical examples[J]. Journal of Biomedical Optics, 2007, 12(4): 041208. doi:  10.1117/1.2772879
    [10] SZKULMOWSKI M, WOJTKOWSKI M. Averaging techniques for OCT imaging[J]. Optics Express, 2013, 21(8): 9757-9773. doi:  10.1364/OE.21.009757
    [11] MARTINEZ-CONDE S, MACKNIK S L, HUBEL D H. The role of Fixational eye movements in visual perception[J]. Nature Reviews Neuroscience, 2004, 5(3): 229-240. doi:  10.1038/nrn1348
    [12] KLEIN T, HUBER R. High-speed OCT light sources and systems [Invited][J]. Biomedical Optics Express, 2017, 8(2): 828-859. doi:  10.1364/BOE.8.000828
    [13] 余轮, 魏丽芳. 眼底图像配准技术研究进展[J]. 生物医学工程学杂志,2011,28(5):1043-1047.

    YU L, WEI L F. Progress of research in retinal image registration[J]. Journal of Biomedical Engineering, 2011, 28(5): 1043-1047. (in Chinese)
    [14] NIEMEIJER M, GARVIN M K, LEE K, et al. Registration of 3D spectral OCT volumes using 3D SIFT feature point matching[J]. Proceedings of SPIE, 2009, 7259: 72591I.
    [15] NIEMEIJER M, LEE K, GARVIN M K, et al. Registration of 3D spectral OCT volumes combining ICP with a graph-based approach[J]. Proceedings of SPIE, 2012, 8314: 83141A.
    [16] CHEN M, LANG A, YING H S, et al. Analysis of macular OCT images using deformable registration[J]. Biomedical Optics Express, 2014, 5(7): 2196-2214. doi:  10.1364/BOE.5.002196
    [17] KHANSARI M M, ZHANG J, QIAO Y CH, et al. Automated deformation-based analysis of 3D optical coherence tomography in diabetic retinopathy[J]. IEEE Transactions on Medical Imaging, 2020, 39(1): 236-245. doi:  10.1109/TMI.2019.2924452
    [18] CHENG J, DUAN L X, WONG D W K, et al.. Speckle reduction in optical coherence tomography by image registration and matrix completion[C]//Proceedings of the 17th International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer, 2014: 162-169.
    [19] LEE S, LEBED E, SARUNIC M V, et al. Exact surface registration of retinal surfaces from 3-D optical coherence tomography images[J]. IEEE Transactions on Biomedical Engineering, 2015, 62(2): 609-617. doi:  10.1109/TBME.2014.2361778
    [20] LEZAMA J, MUKHERJEE D, MCNABB R P, et al. Segmentation guided registration of wide field-of-view retinal optical coherence tomography volumes[J]. Biomedical Optics Express, 2016, 7(12): 4827-4846. doi:  10.1364/BOE.7.004827
    [21] DU X Y, GONG L, SHI F, et al.. Non-rigid registration of retinal OCT images using conditional correlation ratio[M]//CARDOSO M J, ARBEL T, MELBOURNE A, et al.. Fetal, Infant and Ophthalmic Medical Image Analysis. Cham: Springer, 2017: 159-167.
    [22] OJANSIVU V, HEIKKILA J. Image registration using blur-invariant phase correlation[J]. IEEE Signal Processing Letters, 2007, 14(7): 449-452. doi:  10.1109/LSP.2006.891338
    [23] 孙辉, 马天玮. 基于相位相关的目标图像亚像元运动参数估计[J]. 液晶与显示,2011,26(6):858-862. doi:  10.3788/YJYXS20112606.0858

    SUN H, MA T W. Sub-pixel motion estimation based on phase-only correlation[J]. Chinese Journal of Liquid Crystals and Displays, 2011, 26(6): 858-862. (in Chinese) doi:  10.3788/YJYXS20112606.0858
    [24] 万钇良, 王建立, 张楠, 等. 一种基于相位相关与子图像的偏振图像配准方法[J]. 液晶与显示,2019,34(5):530-536. doi:  10.3788/YJYXS20193405.0530

    WAN Y L, WANG J L, ZHANG N, et al. Polarized image registration method based on phase correlation and sub-graph[J]. Chinese Journal of Liquid Crystals and Displays, 2019, 34(5): 530-536. (in Chinese) doi:  10.3788/YJYXS20193405.0530
    [25] PERONA P, MALIK J. Scale-space and edge detection using anisotropic diffusion[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990, 12(7): 629-639. doi:  10.1109/34.56205
    [26] 王亚强, 陈波. 一种改进的各向异性扩散超声图像去噪算法[J]. 液晶与显示,2015,30(2):310-316. doi:  10.3788/YJYXS20153002.0310

    WANG Y Q, CHEN B. Improved anisotropic diffusion ultrasound image Denoising algorithm[J]. Chinese Journal of Liquid Crystals and Displays, 2015, 30(2): 310-316. (in Chinese) doi:  10.3788/YJYXS20153002.0310
    [27] PUVANATHASAN P, BIZHEVA K. Interval type-II fuzzy anisotropic diffusion algorithm for speckle noise reduction in optical coherence tomography images[J]. Optics Express, 2009, 17(2): 733-746. doi:  10.1364/OE.17.000733
  • [1] 周文舟, 范晨, 胡小平, 何晓峰, 张礼廉.  多尺度奇异值分解的偏振图像融合去雾算法与实验 . 中国光学, doi: 10.37188/CO.2020-0099
    [2] 程亚亚, 于化东, 于占江, 许金凯, 张向辉.  微铣刀同轴全息图像增强方法 . 中国光学, doi: 10.37188/CO.2019-0217
    [3] 黄鹤, 李昕芮, 宋京, 王会峰, 茹锋, 盛广峰.  多尺度窗口的自适应透射率修复交通图像去雾方法 . 中国光学, doi: 10.3788/CO.20191206.1311
    [4] 蔡怀宇, 张玮茜, 陈晓冬, 刘珊珊, 韩晓艳.  眼科光学相干层析成像的图像处理方法 . 中国光学, doi: 10.3788/CO.20191204.0731
    [5] 王浩, 张叶, 沈宏海, 张景忠.  图像增强算法综述 . 中国光学, doi: 10.3788/CO.20171004.0438
    [6] 郝建坤, 黄玮, 刘军, 何阳.  空间变化PSF非盲去卷积图像复原法综述 . 中国光学, doi: 10.3788/CO.20160901.0041
    [7] 许廷发, 苏畅, 罗璇, 卞紫阳.  基于梯度和小波变换的水下距离选通图像去噪 . 中国光学, doi: 10.3788/CO.20160903.0301
    [8] 陶小平.  大口径反射镜加工机床在线检测高精度对准方法 . 中国光学, doi: 10.3788/CO.20150806.1027
    [9] 许佳佳.  结合Harris与SIFT算子的图像快速配准算法 . 中国光学, doi: 10.3788/CO.20150804.0574
    [10] 严洁, 阮友田, 薛珮瑶.  主被动光学图像融合技术研究 . 中国光学, doi: 10.3788/CO.20150803.0378
    [11] 王新华, 黄玮, 欧阳继红.  多探测器拼接成像系统实时图像配准 . 中国光学, doi: 10.3788/CO.20150802.0211
    [12] 史光辉, 杨威.  用于图像拼接的电视摄像光学系统 . 中国光学, doi: 10.3788/CO.20140704.0638
    [13] 郑贤良, 刘瑞雪, 夏明亮, 曹召良, 宣丽.  液晶自适应光学视网膜校正成像技术研究 . 中国光学, doi: 10.3788/CO.20140701.098
    [14] 阎春生, 廖延彪, 田芊.  层析成像图像重建算法综述 . 中国光学, doi: 10.3788/CO.20130605.0617
    [15] 杨文波, 马天玮, 刘剑.  非局部变分修复法去除高密度椒盐噪声 . 中国光学, doi: 10.3788/CO.20130606.876
    [16] 孙辉, 李志强.  基于相位相关的匀速直线运动 模糊图像位移参数估计 . 中国光学, doi: 10.3788/CO.20120502.0174
    [17] 于前洋.  “视频图像处理专集”导读 . 中国光学,
    [18] 韩广良.  高频信息矢量匹配实现异源图像配准 . 中国光学,
    [19] 孙辉, 李志强, 孙丽娜, 郎小龙.  一种空域和频域相结合的运动图像亚像素配准技术 . 中国光学,
    [20] 孙辉, 李志强, 孙丽娜, 郎小龙.  基于相位相关的亚像素配准技术及其在电子稳像中的应用 . 中国光学,
  • 加载中
图(8) / 表 (1)
计量
  • 文章访问数:  2
  • HTML全文浏览量:  0
  • PDF下载量:  0
  • 被引次数: 0
出版历程

扫频光学相干层析视网膜图像配准去噪算法

doi: 10.37188/CO.2020-0130
    基金项目:  国家重点研发计划(No.2017YFC0109901);天津市自然科学基金项目(No.15JCQNJC14200)
    作者简介:

    蔡怀宇(1965—),女,湖南涟源人,博士,教授,硕士生导师,1991年、2000年于天津大学分别获得硕士、博士学位,主要从事信息光学、光电技术及仪器和图像处理等方面的研究。E-mail:hycai@tju.edu.cn

    韩晓艳(1996—),女,内蒙古乌兰察布人,天津大学精密仪器与光电工程技术学院硕士研究生,2018年于南京理工大学获得学士学位,主要从事光学相干层析成像方面的研究。E-mail:xyhan@tju.edu.cn

  • 中图分类号: TN247

摘要: 多帧叠加平均处理是去除扫频光学相干层析系统的散斑噪声、获得较为清晰结构信息的有效方法,但眼睛的震颤、漂移、微眼跳等生理特性和系统光路特性会使图像之间存在几何变换,导致叠加效果不佳、结构稳定性差,为此本文提出一种基于灰度分布信息和目标几何信息相结合的配准算法。该方法根据图像平均灰度分布提取包含目标信息的感兴趣区域,通过相位相关算法和基于分段拟合的灰度投影算法的双重作用校正图像的平移变换;通过拟合视网膜上边界作为特征点迭代确定最佳旋转参数,并再次重新估计平移参数,实现图像的刚性配准;最后通过轴向扫描一对一映射法以能量函数为约束条件实现图像的非刚性配准。对活体兔眼进行实验,结果表明,本文算法配准后的叠加图像边界清晰,结构信息增强,信噪比和对比度平均有效提高一倍多。本算法适用于强噪声视网膜B-Scans图像的配准,能满足多种类型OCT系统的叠加成像需要,具有较高的鲁棒性和图像配准精度。

English Abstract

蔡怀宇, 韩晓艳, 娄世良, 汪毅, 陈文光, 陈晓冬. 扫频光学相干层析视网膜图像配准去噪算法[J]. 中国光学. doi: 10.37188/CO.2020-0130
引用本文: 蔡怀宇, 韩晓艳, 娄世良, 汪毅, 陈文光, 陈晓冬. 扫频光学相干层析视网膜图像配准去噪算法[J]. 中国光学. doi: 10.37188/CO.2020-0130
CAI Huai-yu, HAN Xiao-yan, LOU Shi-liang, WANG Yi, CHEN Wen-guang, CHEN Xiao-dong. Speckle noise reduction in swept-source optical coherence tomography by retinal image registration[J]. Chinese Optics. doi: 10.37188/CO.2020-0130
Citation: CAI Huai-yu, HAN Xiao-yan, LOU Shi-liang, WANG Yi, CHEN Wen-guang, CHEN Xiao-dong. Speckle noise reduction in swept-source optical coherence tomography by retinal image registration[J]. Chinese Optics. doi: 10.37188/CO.2020-0130
    • 扫频光学相干层析(Swept Source Optical Coherence Tomography, SS-OCT)是一种高分辨率的非侵入式医学成像技术[1-4],广泛应用于探测视网膜的微观结构,为黄斑病变、糖尿病性视网膜病变、年龄相关性黄斑变性等眼底疾病提供了十分有效的检查方式[5-6]。由于低相干干涉的成像机理,SS-OCT系统会不可避免地引入呈现随机分布的散斑噪声[7],使B-Scan图像信噪比低,结构信息模糊、部分甚至完全缺失,影响了对眼底疾病诊断的准确性。因此,实现在强噪声背景下有效去除散斑噪声、获得结构信息具有重要的研究价值。

      目前,研究者们已提出了很多减少散斑噪声的方法,主要分为单帧去噪方法[8]和多帧叠加平均去噪方法[9]两大类。基于单帧的去噪方法容易平滑图像细节,造成图像信息丢失,无法满足医学图像处理的基本要求。基于多帧的叠加平均去噪方法利用SS-OCT系统在短时间内获取的多帧图像进行叠加平均处理,可以在保证不改变视网膜拓扑结构的基础上有效减少散斑噪声,增强结构信息,已成为去除散斑噪声的最优方法[10]。然而,眼睛的震颤、漂移、微眼跳等生理特性[11]和系统光路的振镜连续扫描模式会导致图像间存在微小的几何变换,直接叠加效果较差且易产生重影区域,最终影响视网膜结构信息的稳定性。因此,在进行叠加平均处理之前必须通过配准算法校准图像的几何变换,从而保证叠加的稳定性要求。图像配准不依赖于高、超高速扫频光源[12]等硬件的支持,可以在低成本的情况下达到高精度校正,并且有助于推动视网膜三维重建、视网膜三维血流检测等领域的发展。

      现有的视网膜配准方法主要分为基于特征的配准方法和基于灰度的配准方法两大类[13]。基于特征的配准方法通过提取图像对的相应特征并计算相应特征的匹配度来校正图像的几何变换。Niemeijer等人[14]提出基于三维尺度不变特征变换特征点的刚性配准方法,得到了效果较好的二维视网膜图像的配准结果。Niemeijer等人[15]又通过提取眼底血管图像的分支和交叉点等特征,利用基于最近点迭代算法和图论的方法实现B-Scans图像的配准。Chen[16]和Khansari[17]等人分别提出了利用黄斑中心凹的特征信息初始配准正常黄斑图像和异常黄斑图像的三维非刚性配准方法。这些方法计算速度快,能够得到较好的匹配结果,但配准精度受限于特征的选择及其匹配精度,易受噪声影响出现错误匹配或无法匹配的问题,特别是在强噪声图像中易出现配准精度不高、配准效果不理想的现象。基于灰度的配准方法利用图像的灰度信息并采用不同的相似度度量标准求解图像间的最佳变换关系。Chen等人[18]提出使用全局对齐和基于快速迭代菱形搜索的局部对齐来配准图像,但易受强噪声干扰。Lee等人[19]提出基于视网膜分层信息的层匹配算法,以基于灰度信息的分割算法生成的视网膜层表面模型为标准实现宽视场图像[20]的配准,但配准精度很大程度取决于分割模型的精确性,尤其不适用强噪声背景下视网膜结构层模糊、缺失的图像。Du等人[21]提出以相邻图像空间分布的条件相关比信息代替互信息作为相似性度量的非刚性配准方法,但当图像的噪声、结构等因素有较大变化时,存在配准结果不精确的问题。综上可知,现有的视网膜配准算法对特征信息(不变特征点、眼底血管图像的分支交叉点或黄斑中心凹等特殊区域等)或灰度信息(视网膜分层信息等)具有较强的依赖性,容易出现配准信息无法运用或配准不精确的问题,不能有效校正多种类型OCT设备采集的视网膜图像的几何变换。

      针对现有算法对强噪声图像配准不精确的问题,本文根据眼睛生理特性和系统光路特性分析了图像的几何变换,并提出了一种基于灰度分布信息和目标几何信息相结合的高精度配准方法,该算法不过分依赖于图像的特征信息,实现了对强噪声背景下目标信号的精确配准;对配准后的图像进行叠加平均处理,能够在稳定视网膜结构信息的同时有效减少散斑噪声。与其他算法相比,本文算法的鲁棒性较高,对受不同程度噪声污染的视网膜图像都能有效配准,且配准精度较高,能够保证图像叠加平均处理的有效性,大幅度提高视网膜图像的质量;也可以减少硬件成本,有助于推动智能医疗和临床医学等领域的发展。

    • 在SS-OCT系统采集视网膜的过程中,通过视线引导装置进行眼睛的定位定焦,可以快速、准确地对眼底区域进行高分辨率的成像。但在注视的过程中,眼睛会有漂移、震颤、微眼跳三种微弱的生理活动[11]。漂移是一种沿视轴移动的缓慢运动,具有较小的振幅范围;震颤是一种在最小振幅范围内无规则的高频运动,具有平移和旋转的双重特性;微眼跳是一种沿垂直视轴方向的跳跃式运动,具有最大的振幅范围和较小的振动频率。因此,眼睛的运动可以归结为沿着视轴方向、垂直视轴方向的平移和多方向旋转运动,如图1(a−c)所示。这些运动会不同程度地影响系统对视网膜区域的采集,降低多帧图像直接叠加去噪的有效性。

      图  1  OCT序列图像几何变换原因分析(a)眼睛轴向运动(b)眼睛垂轴运动(c)眼睛旋转运动(d)多角度光线入射

      Figure 1.  Analysis of reasons for geometric transformation of OCT sequence images (a) Axial eye movement (b) Vertical eye movement (c) Rotational eye movement (d) Multi-angle light incidence

      假设测量时,SS-OCT系统的光轴与人眼的视轴完全重合,每帧图像获取时间极短。在获取系列图像过程中,眼睛沿着轴向方向移动(图1(a)),则会引起轴向扫描起始点的沿轴移动,在视网膜各层间隔保持不变的情况下造成图像的整体轴向平移;当眼睛沿着垂直视轴方向移动时(图1(b)),系统对视网膜的采样起始点也会随之上下移动,造成图像的整体垂轴平移。因此,眼睛的轴向、垂轴运动使图像之间存在平移变换,需要通过平移参数估计实现对图像间的重叠区域的配准叠加处理。眼睛的旋转现象是一种复杂的生理活动,其运动过程不具规则性。为了简化研究,不考虑眼睛旋转时(图1(c))引起的成像区域错动,仅对图像序列中视网膜层状结构的微小旋转角度进行处理,通过旋转参数估计实现图像的去旋转变换。

      OCT系统的光路特性也会引起图像的畸变。SS-OCT系统采用振镜连续扫描模式,将不同角度的平行光聚焦于眼底不同区域,通过连续采集得到一帧完整的二维视网膜B-Scan图像。当光线沿光轴垂直入射眼睛时,光线方向为视网膜层的法线方向,实测视网膜层间距具有最短距离;当光线以一定角度入射眼睛时,光线方向偏离视网膜层法线方向,实测层不是层间最短距离,会受入射角度影响呈现非线性变换关系,即造成图像中每行扫描线产生了畸变。扫描角度越大,畸变效应越严重,需要通过非线性估计进行校正。

      综上所述,眼睛的漂移、震颤、微眼跳等生理现象和系统光路特性会影响SS-OCT系统对视网膜的成像,造成图像之间存在平移、旋转、畸变等几何变换,在叠加平均处理之前必须对多帧图像进行针对性配准处理。

    • SS-OCT系统获取的原始单帧视网膜B-Scan图像包含500条A-Scan扫描线,每行为1条A-Scan轴向扫描线。原始图像存在大量的散斑噪声,大部分是不具有所需信息的背景区域,视网膜结构信息较弱,只具有相对较明显的视网膜上边界。根据图像特征和图像几何变换分析,本文配准算法主要分为预处理、刚性配准和非刚性配准三部分。首先,根据图像的列平均灰度分布划分感兴趣区域和背景区域,以减少后续配准运算时间。然后,利用基于灰度信息的相位相关(Phase Only-Correlation, POC)算法校正图像的平移变换,并通过计算图像间的行灰度投影曲线的匹配度进一步求解垂轴平移量,实现图像的平移配准;通过拟合视网膜的上边界作为图像的几何特征点迭代匹配确定最佳旋转参数,并再次按平移方法搜索最佳平移量,实现图像的刚性配准。最后,以参考图像的上边界为标准在最小能量函数的约束条件下对目标图像进行轴向微平移,校正边缘区域的畸变,实现图像的非刚性配准。算法流程如图2所示。

      图  2  配准算法流程图

      Figure 2.  Flow chart of the registration algorithm

    • SS-OCT系统采集的视网膜B-Scan图像包括以视网膜信号为主体的感兴趣区域(Region of Interest, ROI)和大部分无效暗背景区域(Region of Background, ROB),背景区域过大会降低图像配准的速度。因此,本文利用图像的列平均灰度分布规律提取感兴趣区域,即先计算各列像素的平均灰度,再以图像的平均灰度值为阈值划分感兴趣区域和背景区域;为保证上边界的完整性,采用降低左阈值法对ROI区域进行微调整。预处理过程如图3所示,图3(b)为原始图像(图3(a))的列平均灰度分布,确定阈值后得到视网膜感兴趣区域(图3(c))。考虑成像系统和图像几何变换因素,在图3(c)中,画出扫描角度小于$ \pm {4.5}^{ \circ }$$ {R}_{H}$区域为视网膜结构信息较为稳定的中央区域,超出部分则为视网膜信号较弱的边缘区域($ {R}_{L}$区域)。

      图  3  预处理过程(a)原始视网膜B-Scan图像(b)列平均灰度分布图(c)ROI图像

      Figure 3.  Pretreatment process (a) Original B-Scan retinal image (b) Column average gray distribution map (c) ROI image

    • 眼睛的轴向、垂轴运动引起轴向扫描的全局变化,造成视网膜图像之间存在平移变换。由于图像具有散斑噪声多、干涉信号弱等明显的特点,导致大部分特征点配准算法难以应用,因此本文采用抗噪声较好的基于灰度信息的POC算法对图像进行轴向配准,可以依靠视网膜上边界的高灰度信息准确估计轴向平移量$ \Delta y$;针对图像在垂轴方向纹理信息较弱的问题,提出基于分段拟合的行灰度投影算法进行校正,通过与POC算法的双重作用精确估计垂轴平移量$ \Delta x$

      POC算法[22-24]利用基于傅里叶变换的平移不变性将图像在空域中的平移量转化成对POC函数峰值位置的求取。假设目标图像$ {f}_{2}\left(x,y\right)$是由参考图像$ {f}_{\rm{1}}\left(x,y\right)$平移得到,对二者进行离散傅里叶变换得到$ {F}_{\rm{1}}\left(\rm{u},v\right)$$ {F}_{\rm{2}}\left(\rm{u},v\right)$,计算二者互功率谱的相位$ G\left(u,v\right)$,再进行反傅里叶变换可得到POC函数:

      $$ p\left(x,y\right)\approx \frac{\mathrm{sin}\left[\pi \left(x-\Delta x\right)\right]}{\pi \left(x-\Delta x\right)}\frac{\mathrm{sin}\left[\pi \left(y-\Delta y\right)\right]}{\pi \left(y-\Delta y\right)}$$ (1)

      其中,$ \left(x,y\right)$为空域图像的像素坐标,$ \left(u,v\right)$是频域图像的像素坐标,$ p\left(x,y\right)$为POC函数的最大峰值。通过计算$ p\left(x,y\right)$峰值的位置就可确定两图像间的空域平移量$ \left(\Delta x,\Delta y\right)$

      针对图像垂轴信号较弱的问题,对于POC配准后的多帧图像,采用行灰度投影算法进行二次配准,可以高速、高精度地估计垂轴平移量,避免了对图像纹理信息和边界形状的依赖性。行灰度投影算法是对图像的行像素进行投影,得到灰度映射曲线,相应的行灰度投影值$ m\left(x\right)$的计算公式为:

      $$ m\left(x\right)=\frac{1}{C}{\displaystyle \sum _{y=1}^{C}I\left(x,y\right)}$$ (2)

      其中,$ C$为最大列数,$ I\left(x,y\right)$是坐标$ \left(x,y\right)$处的像素灰度值。当目标曲线和参考曲线的峰点横坐标达到最大程度的重合,即曲线对的重叠面积达到最大时,目标图像才能实现最优垂轴配准。为精确估计垂轴平移量,先利用参考曲线谷点集分割曲线对,设定曲线段长度均值为阈值进行合并处理;通过旋转位于曲线上方的包络线重新拟合曲线段,使多峰高度一致;再对曲线段做归一化处理即可。相应的曲线分段、拟合过程如图4所示,处理后的曲线对消除了纵坐标差值的影响,可以通过精确计算曲线对的重叠面积估计垂轴平移量。

      图  4  分段拟合处理过程(a)曲线分割结果(黑色虚线为分割线)(b)原始曲线段(最佳连接点A、黑色直线为包络线)(c)处理后的曲线段

      Figure 4.  Fitting process of curve segments (a) Results of curve segmentation (the dotted black line is the segmentation line) (b) Original curve segment (optimal connection point A, the straight black line is the envelope) (c) Curve segment after processing

    • 针对眼睛绕轴转动导致图像间产生微小旋转变换的问题,本文提出基于视网膜上边界的特征点迭代方法。视网膜上边界是具有高灰度值的几何信息,通过局部搜索策略进行精确提取;以目标边界作为旋转中心特征点集合,采用迭代方式确定图像的最佳旋转参数;按平移配准方法再次确定最佳平移量,实现图像的刚性配准。

      为精确提取强噪声图像中的视网膜上边界,先采用各向异性扩散滤波算法[25-26]对图像进行滤波,可以在保留边界细节的同时降低噪声干扰。在边缘区域内,边界信号明显强于视网膜内部有效信号和背景信号,通过局部搜索最大灰度像素能够有效确定上边界一点,并以该点为基础采用局部行搜索范围初步确定上边界集合$ \left\{\left({x}_{m},{y}_{m}\right)\right\}$。计算相邻特征点的列坐标差平均值:

      $$ md=\frac{{\displaystyle \sum _{m=1}^{n\rm{-1}}\left|{y}_{m+\rm{1}}-{y}_{m}\right|}}{n-1}$$ (3)

      其中$ n$为总点数,根据统计学规律可知,当满足$ \left|{y}_{m+\rm{1}}-{y}_{m}\right|<0.95\rm{md}$条件时剔除错误点,并采用双线性插值法得到完整、准确的特征点集合。图5(a)是上边界提取过程的示意图,A点为下边缘区域的最大灰度像素点,红色单箭头线为行局部搜索范围,红色曲线为拟合后的视网膜上边界。

      平移配准后的图像(图5(b))在中央区域具有较好的配准效果,上边界吻合度较高,通过POC算法和灰度投影算法的双重作用有效校正了该区域的平移变换。因此,以目标中央区域内的边界特征点为旋转中心集合,设定初始旋转角度为$ {\theta }_{\rm{0}}$,旋转角度$ \theta $以一定间隔$ \mu $作为步长因子不断改变$ \left(\theta =\mu {\theta }_{\rm{0}}\right)$,采用互信息的评价方式求解不同旋转参数下图像对的相似性,确定目标图像的最优参数,实现图像的去旋转变换;最后通过迭代平移配准方法再次估计最佳平移量。相较于平移配准后的叠加图(图5(b)),旋转配准方法明显提高了叠加图上边界的一致性(图5(c))。

      图  5  刚性配准结果(a)上边界提取过程(b)平移配准后的叠加图(c)旋转配准后的叠加图(绿色为参考图像,紫色为目标图像)

      Figure 5.  Results of rigid registration (a) Upper boundary extraction process (b) Averaging image after translation registration (c) Averaging image after rotation registration (Green is the reference image, purple is the target image)

    • SS-OCT系统发出的近轴光线与光轴的夹角十分微小,在一定的角度范围内可认为是垂直入射到角膜顶点处,可以忽略其对轴向扫描光程的影响,能够通过刚性配准有效校正中央区域的几何变换;但随着扫描角度的增大,轴向扫描光程在系统光学结构的影响下会出现不同程度的变化,导致图像中视网膜各层间距发生变化。因此,需针对图像的边缘区域进行局部畸变校正,方法是采用图像对的轴向扫描一对一映射进行非刚性配准变换。

      进行配准时,只考虑轴向平移,避免因径向平移导致轴向扫描的重叠和视网膜拓扑结构的改变。为了进一步提高配准的精度,先利用一维高斯滤波器对轴向扫描进行平滑滤波,减少散斑噪声的影响;在轴向扫描一对一映射的过程中,设置能量函数$ E\left(\Delta y\right)$作为约束条件,即:

      $$ \left\{ \begin{array}{l} \min E\left( {\Delta y} \right) = \displaystyle\sum\limits_{y = 1}^C {\left| {{I_o}\left( {x,y + \Delta y} \right) - {I_r}\left( {x,y} \right)} \right|} \\ - 2D {\text{≤}} \Delta y {\text{≤}} 2D \end{array} \right.$$ (4)

      其中,$ {I}_{o}\left(x,y+\Delta y\right)$是目标图像在坐标$ \left(x,y+\Delta y\right)$处的灰度值,$ {I}_{r}\left(x,y\right)$是参考图像在坐标$ \left(x,y\right)$处的灰度值,$ D$为图像对上边界之间的相对距离。能量函数反映得是一对一轴向扫描的灰度分布差异,在一定的轴向平移范围$\left(-2D{\text{≤}} \Delta y{\text{≤}} 2D\right)$内寻找$ E\left(\Delta y\right)$的最小值,得到高精度的轴向匹配结果。非刚性配准如图6(b)所示,相较于刚性配准后的叠加图(6(a)),非刚性配准校正了图像的边缘畸变,避免了重影区域的产生,提高了图像的配准精度。

      图  6  非刚性配准结果(a)刚性配准后的叠加图及局部放大图(b)非刚性配准后的叠加图及局部放大图(绿色为参考图像,紫色为目标图像)

      Figure 6.  Results of non-rigid registration (a) Average image obtained after rigid registration (b) Average image obtained after non-rigid registration (The red dashed box represents an enlarged partial view; Green is the reference image, and purple is the target image)

    • 实验采用自主搭建的视网膜成像SS-OCT系统进行,系统光源中心波长为1060 nm、带宽为70 nm、频率为10 kHz(Santec HSL-1),在组织内的实际纵向分辨率为9.53 um,实际成像深度为10.83 mm。对活体健康兔眼眼底进行成像(活体兔为新西兰品种,年龄4岁,体重2.5 kg),实验过程遵守动物医学伦理要求,得到相关部门审批。通过对SS-OCT系统采集到的多帧B-Scan图像进行配准、叠加处理,得到如图7所示的结果。相较于原始单帧图像(图7(a)),配准后的叠加图像(图7(b))的有效信号明显增强,散斑噪声得到有效去除,且具有清晰、一致的视网膜上边界。

      图  7  本算法实验结果图(a)原始单帧图(b)配准后的叠加图

      Figure 7.  The experimental results of the algorithm (a) Original single frame image (b) Averaging image after registration

      为了评价本文方法的性能,采用信噪比(signal-to-noise ratio, SNR)和对比度(contrast-to-noise ratio, CNR)的客观标准来定量分析图像质量,相应的计算公式[27]分别为:

      $$ SNR=\rm{10}{\mathrm{log}}_{10}\left(\frac{\mathrm{max}\left[I{\left(x,y\right)}^{2}\right]}{{\sigma }_{b}^{2}}\right)$$ (5)
      $$ CNR=\frac{1}{R}\left(\frac{{\displaystyle \sum _{k=1}^{K}\left({\mu }_{e}-{\mu }_{b}\right)}}{\sqrt{{\sigma }_{e}^{2}+{\sigma }_{b}^{2}}}\right)$$ (6)

      式中,$ {\mu }_{b}$$ {\sigma }_{b}$分别为背景区域的平均值和标准差,$ {\mu }_{e}$$ {\sigma }_{e}$分别为第$ k$个视网膜有效区域的平均值和标准差,$ K$为有效区域的总个数。以图8(a)所示的原始单帧图像为例,选定一个较大的黄色矩形框为背景区域和七个红色矩形框为视网膜有效区域($ K=10$)。利用本文方法、文献[18]方法和文献[21]方法对多帧连续图像进行配准,计算不同帧数叠加图在规定区域下的SNR和CNR值,相应的变化趋势分别如图8(b)图8(c)所示。

      图  8  SNR和CNR随叠加帧数递增的曲线变化图(a)区域选择(b)SNR(c)CNR

      Figure 8.  SNR and CNR curve change graph with an increasing number of multiple frames (a) Area selection (b) SNR (c) CNR

      图8(b)8(c)可以看出,随着叠加帧数的递增,本文方法、文献[18]方法和文献[21]方法的SNR和CNR都不断提高,但本文方法的结果始终高于其他两种方法的结果,相较于文献[18]方法,SNR平均提高了1.12倍,CNR平均提高了1.19倍;相较于文献[21]方法,SNR平均提高了1.50倍,CNR平均提高了1.48倍,如表1所示。因此,本文方法是一种高精度的图像配准方法,可以保证对强噪声图像的配准精度,从而能通过多帧叠加平均处理在增强视网膜结构信息的同时有效去除散斑噪声。

      表 1  配准算法对比结果

      Table 1.  Comparison of the registration algorithms

      比较对象SNR提高倍数CNR提高倍数
      本文方法/文献[18]方法1.121.19
      本文方法/文献[21]方法1.501.48
    • 通过分析眼睛的生理特征和系统的光路特性可知图像存在平移、旋转、畸变等几何变换;针对不同的几何变换问题,提出了一种基于灰度分布信息和目标几何信息相结合的强噪声视网膜图像配准算法,包括目标区域提取的预处理、全局的刚性配准和局部的非刚性配准三部分。将该方法应用在自主搭建SS-OCT系统采集的视网膜B-Scan图像的配准中,配准后的叠加图像的成像质量明显得到提高,结构信息显著增强,边界吻合度具有较高的一致性,且与已有算法相比,性能指标(信噪比和对比度)提升一倍多。因此,本文算法不过分依赖于图像的特征信息,对不同程度噪声污染的视网膜图像具有高精度的配准效果,通过叠加平均处理可以在减少硬件成本的同时有效提高眼底疾病的诊断准确率,可以促进眼科医学领域的发展,并且也有助于视网膜三维重建、视网膜三维血流检测等领域的发展。

参考文献 (27)

目录

    /

    返回文章
    返回