OCT图像运动伪差校正算法的研究

【作者】张明蓉，李跃杰,2，王立伟，高振玉

1 北京协和医学院 中国医学科学院生物医学工程研究所，天津市，300192

2 天津市眼科医学设备技术工程中心，天津市，300384

OCT图像运动伪差校正算法的研究

【作者】张明蓉1，李跃杰1,2，王立伟1，高振玉1

1 北京协和医学院 中国医学科学院生物医学工程研究所，天津市，300192

2 天津市眼科医学设备技术工程中心，天津市，300384

目的 光学相干层析成像系统（Optical Coherence Tomography, OCT）进行动态成像时，运动伪差会使图像产生畸变，从而可能导致临床诊断的误诊以及错诊，运动伪差的校正已经是亟待解决问题。方法 首先用改进的自适应非线性复合扩散滤波器对连续采集的序列图像进行预处理降噪，然后用图像的质量重心对畸变的体数据在Y方向进行粗略的对齐，再用差值平均法对Y方向的运动伪差进行校正。结果 经过校正，纵切2D图像毛刺消失，表面平滑；3D图像表面平滑，结构清晰，层次分明，视网膜en-face单层切片图像清晰，可以观察到眼底组织细微结构。结论 经过运动伪差校正后，眼底视网膜的生理结构得到了真实的反应。

光学相干层析成像；运动伪差；校正

0 引言

光学相干层析成像(Optical Coherence Tomography，OCT)技术可以对组织的内部结构和形态进行二维和三维可视化[1]。具有高分辨率、高灵敏度、活体、实时、非侵入、横向分辨率和纵向分辨率相互独立等优点。自从OCT出现，它已经成了眼科疾病诊断和研究的一种标准工具[2]。频域OCT的出现[3-4]极大的提高了成像的灵敏度和扫描速度[5-10]，使OCT系统实时成像成为可能。

但是，OCT系统（甚至一些超高速的OCT系统）在获取眼底图像时，由于受试者眼睛的运动（包括人眼的飘逸、眨眼、扫视以及头部等不自主的运动）会引起运动伪差。运动伪差会使图像发生畸变，图像的畸变会使数据丧失连续性，这样三维重建出的三维图像将会丢失一些组织的结构信息，使得视网膜的生理结构不能得到真实的反应。实时成像过程中由于OCT系统的高度敏感性，甚至会导致成像的失败，从而导致一些眼科疾病的误诊。所以图像运动伪差的校正在OCT技术中占据相当重要的地位，已经成为了相关领域的研究热点。进年来一些研究小组做了相关的研究，提出采用硬件方法实现运动伪差校正，其中应用硬件来消除运动伪差最典型的方法为2012年Vienola KV[11]发表的结合追踪激光扫描检眼镜法 (TSLO) 进行补偿来追踪实时的眼运动。目前一些商业OCT系统采用硬件辅助实现眼睛运动的实时校正，如德国海德堡公司SPECTRALIS OCT眼科OCT系统，美国加州Optovue公司的RTVue OCT系统，Physical Sciences公司的追踪OCT系统[12-14]。同时也有相关研究采用软件算法实现运动伪差校正，但是都只能解决局部的运动伪差[10, 15-20]。

本文在不增加任何硬件的前提下，设计伪差校正算法对数据采集过程中产生的运动伪差进行校正，并且得到了很好的结果。

1 系统组成

系统结构图如图1所示，采用课题组研制的3D眼科前后节FD-OCT光学系统进行光学检测，通过控制信号扫描方式采集信号，在采集信号时可能由于一些不可控因素导致产生运动伪差，所以对信号进行运动伪差处理。首先进行预处理，对原始数据用改进的非线性复合扩散滤波器进行预处理去噪；然后结合慢轴方向快扫的图像应用差值平均法对体数据Y方向进行校正；最后进行3D可视化显示。

图1 系统结构图Fig.1 The system architecture

2 运动伪差的校正

2.1 信号扫描

OCT系统体数据扫描模式如图2(左)所示，图2(右)是3D数据的虚拟切割平面。在这个标准的体数据模型中快轴和慢轴扫描分别对应水平和垂直扫描仪。快轴每扫描一次即为一个B-scan，将获取的原始数据经过后处理以及畸变校正算法之后，进行三维重建，将重建之后的体数据可以虚拟的切出慢轴方向的B-scan和C-scan。

图2 OCT系统体数据扫描模式Fig.2 The OCT scanning pattern in volumetric acquisition

但是在OCT进行动态活检时，由于一些人眼的飘逸、眨眼、扫视以及头部等不自主的运动极易引起运动伪差，其原理如图3所示。

一般可以通过提高图像的扫描速度来降低运动伪差，而运动伪差在很大程度上取决于受试者，在采集数据时，随着扫描时间的增加运动伪差也会增大。在理想状态下，即体数据获取速度（小于100 ms）较快时会达到无运动伪差的模式。但是这种情况在目前的眼科OCT设备是不可能实现，除非减少体数据中A扫描的数量。而减少A扫描的数量，则会严重的限制眼底图像的细节，3D视网膜数据的获取和横向的抽样密度之间相互制约。在扫描时可以像其他眼科成像设备一样固定下巴和前额减少头部运动，从而降低轴向的运动伪差。使受试者凝视一个固定目标也可降低成像过程中的运动伪差。但是，这些方法都不能完全的消除运动伪差，所以在OCT成像过程中运动伪差的校正是非常必须的。

图3 眼睛和头运动产生畸变的原理图Fig.3 Schematic of the eye and head motion caused artifact

2.2 运动伪差的校正算法

首先对获取的原始体数据用改进的非线性复合扩散滤波器进行预处理降噪，然后用图像的质量重心对畸变的体数据在Y方向进行粗略地对齐，然后结合慢轴方向快扫的25帧图像应用差值平均法对体数据Y方向进行精确地校正，其操作原理图如图4所示。

图4 运动伪差校正原理图Fig.4 Schematic of the Motion artifact correction

2.2.1 数据的预处理

课题组研制的3D眼科前后节FD-OCT系统的可达到每秒钟70000个B-scan的成像速度，单帧图像的成像速度非常快，单帧的B-scan图像几乎不会出现畸变，所以本文只对连续的B-scan进行了校正，没有对单个的B-scan进行校正。

在OCT成像系统中，图像的噪声主要为散斑噪声和背景噪声，造成这些噪声的主要因素包括光源光斑大小、时间相干性、传输光的多次散射、相位的畸变以及目标的光学性能等。为了后续处理得到更好更准确的结果，对原始图像进行降噪处理是非常重要的。本文采用Rui Bernardes于2010年提出的改进自适应非线性复合扩散滤波器，对OCT视网膜图像具有更好的噪声消除和层状结构边缘信息保持能力。

非线性各向异性复合扩散是基于式(1)的求解：

扩散系数D定义为：

其中k为阈值参数，扩散系数D通过调节参数k可画出一组曲线，如图5所示。

图5 扩散系数D的曲线Fig.5 Curves of diffusion coefficient D

由图5可看出随着k的变化，扩散系数D在其最大值处的延展范围也发生了变化。随着参数ΔI由低强度区域（较高的k值）到高强度区域（较低的k值）的变化，D区别变得更加明显。由于眼底图像的背景区域属于低强度区域，视网膜区域属于高强度区域。结合扩散系数D的曲线变化规律，改进参数k，促进低强度区域的扩散，降低高强度区域的扩散，即可很好的保留图像的边缘信息。

一般非线性复合扩散算法中的时间步长都是设定的固定值，但是图像梯度对扩散系数影响极大，而且在扩散迭代初始时噪声会使梯度值更大，进而会造成扩散系数的变化，所以将时间步长设置为自适应的，自适应时间步长如式(3)所示：

其中a+b≤1。

2.2.2 校正Y方向的伪差

在校正Y方向的畸变错位时，首先求取每帧图像的质量重心，设图像大小为m×n，i为图像的行数，j为图像的列数，则图像质量重心的定义如下：

根据式(4)求得每帧图像的质量重心之后，以第一帧图像为参考帧，使其质量重心Yc与其后一帧图像的质量重心Y'c做差得到Y0，再以此差值Y0将第二帧图像在Y方向做相应的平移对齐，依此类推每帧图像都以前一帧图像为参考图像进行Y方向的粗略对齐，这样做的目的是为了后续处理工作更加精确。

为了校正Y方向的运动伪差，在获取原始的体数据时，首先沿着慢轴扫描方向快速的获取25帧纵向的B-scan，分别包括轴向的0、1/4、2/4、3/4、4/4的5个位置各扫5帧图像，其对应位置如图6所示。

图6 慢轴5个位置的快扫图Fig.6 B-scan along slow axis in the five position

第一步，将每个位置的5帧纵扫图像叠加平均得到更好的参考帧图像。慢轴方向经过快速扫描得到的原始图像经过预处理后如图7(a)所示，将同一位置的5帧图像叠加平均之后得到的图像如图7(b)所示。由图7可以看出经过叠加平均后的图像信息更加清晰，同时也降低了一些背景噪声。这为后续的运动伪差校正奠定了良好的基础。

图7 五帧图像的叠加平均Fig.7 Five-frame mean

第二步，在原始体数据中的上述5个位置纵切得到5帧图像。此时纵切得到的5帧图像与纵扫的图像所对应的信息量是完全一致的，只是在每列像素Y方向上存在位移差。然后寻找每个位置对应的纵切和纵扫叠加平均的图像IS/OS边界（光感受器内节/外节 inner and outer photoreceptor segment, IS/OS）。寻找图像IS/OS边界时，先将图像进行二值化，再去除面积点数在100以下的区域，再将图像进行腐蚀、膨胀，然后将图像的每一列从下向上寻找第一个由1变为0的值的点，即IS/OS初始边界，为了找到更加准确平滑的边界，从上述初始边界第6个点开始，每一个点与其前5个点坐标均值作比较，若差的绝对值大于4，则用此坐标均值代替此点的坐标，最终处理结果为IS/ OS边界。最后将寻找到的纵切和纵扫叠加平均的图像IS/OS边界的值做差值，为了降低单一量的随机误差，对所得到的5组差值求平均值，此平均值即为每帧所对应图像在Y方向的位移量，以此值将每帧图像进行平移，即可使原始数据在Y方向也得到校正。

3 实验结果与讨论

为了验证上述方法的有效性，采用课题组研制的3D眼科前后节FD-OCT系统随机的采取眼底视网膜中央凹数据，在数据采集过程中并没有对受试者进行强迫固定，而是让受试者在一个轻松自然的状态下进行数据的采集。

采集的眼底视网膜区域大小为12 mm×12 mm，原始数据大小为512×256×682，即256帧图像，每帧包括512个A线，深度为682的图像，将原始数据进行可视化得到3D图像。获取的原始数据为B扫描模式，将其进行纵向切割得到慢轴方向的B-scan，将图像进行en-face单层切片可得到C-scan图像，如图8所示，其中(a)为原始数据重建之后得到的3D眼底视网膜图像，(b)为经过校正后重建的3D眼底视网膜图像。

由图8可以看出，在数据采集过程中由于一些眼球不自主的运动产生了位移，使得重建出来的3D图像的表面凹凸不平，使原来光滑的视网膜表面出现严重的褶皱，不能够真实地反应眼底视网膜的生理结构。经过伪差校正之后，重建的3D图像表面光滑，几乎与真实的眼底视网膜图像没有差别。

图9(a)为原始体数据中央凹处（图7，2/4处）进行纵切得到的慢轴方向的B-acan图像；图9(b)校正后体数据中央凹处（图6，2/4处）进行纵切得到的慢轴方向的B-scan图像；图9(c)为原始体数据深度357处的视网膜en-face单层切片图像；图9(d)校正之后提数据357处的视网膜en-face单层切片图像。

由图9可以看出，由于采集数据时的一些不可控因素使得原始数据发生畸变，图9(a)中纵切B-scan图像表面层次不齐，Y方向错位；图9(c)中中央凹表面参差不齐，并包含错误信息，使得视网膜en-face单层切片图像不能得到真实的反应。由图9(b)和图9(d)可以看出，经过差值平均法校正之后，Y方向的位移也得到很好地校正，毛刺消失，表面平滑，与最初在慢轴方向的快扫图像几乎完全一致。视网膜en-face单层切片图像能清晰地观察到视网膜血管、微血管、黄斑等眼底组织细微结构，更加真实的反应了眼底视网膜中央凹附近的生理结构。

4 结论

本文利用标准FD—OCT系统，在没有增加任何的追踪定位系统以及硬件的情况下，单纯的使用软件编程对OCT系统在动态成像过程中运动伪差进行校正。

B-scan能够清晰地显示视网膜各层细微结构及病理改变，主要应用于视网膜疾病和青光眼等疾病诊断。en-face单层切片图像直观地显示眼底视网膜血管和黄斑等组织的结构信息，主要应用于黄斑病变、眼底神经组织结构变化等疾病诊断。运动伪差校正算法使得畸变纵切的B-scan和en-face单层切片图像得到很好的校正，真实再现眼底视网膜图像，为临床眼科疾病的诊断提供依据，提高了临床诊断的准确性和科学性，同时也为实现OCT动态成像中眼运动的实时追踪校正奠定了基础。

[1] Huang D, Swanson EA, Lin CP, et al. Optical coherence tomography[J]. Science, 1991, 254(5035): 1178-1181.

[2] Drexler W, Fujimoto JG. State-of-the-art retinal optical coherence tomography[J]. Prog Retin Eye Res, 2008, 27(1): 45-48.

[3] Fercher AF Hitzenberger CK, Kamp G, et al. Measurement of intraocular distances by backscattering spectral interferometry[J]. Opt Commun,1995, 117(1): 43-48.

[4] Wojtkowski M, SrinivasanVJ, Ko TH, et al. Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation[J]. Opt Express, 2004, 12(11): 2404-2422 .

[5] Leitgeb R, Hitzenberger C. Performance of fourier domain vs. time domain optical coherence tomography[J]. Opt Express, 2003, 11(8): 889-894.

[6] de Boer JF, Cense B, Park BH, et al. Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography[J]. Opt Lett, 2003, 28(21):2067-2069.

[7] Choma M, Sarunic M. Izatt J. Sensitivity advantage of swept source and Fourier domain optical coherence tomography[J]. Opt Express, 2003, 11(18): 2183-2189.

[8] Nassif N, Cense B, Park B. In vivo high-resolution video-rate spectral-domain optical coherence tomography of the human retina and optic nerve[J]. Opt Express, 2004, 12(3): 367-376.

[9] Nassif N, Cense B. In vivo human retinal imaging by ultrahighspeed spectral domain optical coherence tomography[J]. Opt Lett, 2004, 29(5): 480-482.

[10] Potsaid B, Gorczynska I. Ultrahigh speed spectral / Fourier domain OCT ophthalmic imaging at 70,000 to 312,500 axial scans per second[J]. Opt Express, 2008,16(19): 15149-15169.

[11] Vienola KV, Braaf B, Sheehy CK, et al. Real-time eye motion compensation for OCT imaging with tracking SLO[J]. Biomed Opt Express, 2012, 3(11): 2950-2963.

[12] Braaf B, Vienola KV. Real-time eye motion correction inphaseresolved OCT angiography with tracking SLO[J]. Biomed Opt Express, 2012, 3(11): 2950-2963.

[13] Hammer D, Ferguson RD. Advanced scanning methods with tracking optical coherence tomography[J]. Opt Express, 2005, 13(20): 7937-7947.

[14]Vogel CR, Arathorn DW, Roorda A, et al. Retinal motion estimation in adaptive optics scanning laser ophthalmoscopy[J]. Opt Express, 2006, 14(2): 487-497.

[15] Zawadzki RJ, Fullerb AR, Choi SS, et al. Correction of motion artifacts and scanning beam distortions in 3D ophthalmic optical coherence tomography imaging[J]. Proc SPIE 6426, 2007, doi:10.1117/12.701524.

[16] Antony B, Abràmoff MD, Tang L, et al, Automated 3-D method for the correction of axial artifacts in spectral-domain optical coherence tomography images[J]. Biomed Opt Express, 2011, 2(8): 2403-2416.

[17]Kraus MF, Potsaid B, Mayer MA, et al. Motion correction in optical coherence tomography volumes on a per A-scan basis using orthogonal scan patterns[J]. Biomed Opt Express, 2012, 3(6): 1182-1199.

[18]Ricco S, Chen M, Ishikawa H, et al. Correcting motion artifacts in retinal spectral domain optical coherence tomography via image registration[J]. Med Image Comput Comput Assist Interv, 2009, 12(Pt 1): 100-107.

[19]Ortiz S, Siedlecki D, Grulkowski I,et al. Optical distortion correction in Optical Coherence Tomography for quantitative ocular anterior segment by three-dimensional imaging[J]. Opt Express, 2010, 18(3): 2782-2796.

[20]Hammer D, Ferguson RD, Iftimia N, et al. Advanced scanning methods with tracking optical coherence tomography[J]. Opt Express, 2005, 13(20): 7937-7947.

Study of Motion Artifacts Correction Algorithm in Optical Coherence Tomography Images

【Writers】ZHANG Mingrong1, LI Yuejie1,2, WANG Liwei1, GAO Zhenyu1

1 Peking Union Medical College, Chinese Academy of Medical Sciences Institute of Biomedical Engineering, Tianjin, 300192
2 Tianjin Ophthalmic Medical Device Technology Engineering Center, Tianjin, 300384

Objective Optical coherence tomography images may be distorted by motion artifacts in dynamic imaging, so it may lead to misdiagnosis in clinical diagnosis. Motion artifacts correction has become an urgent issue in optical coherence tomography imaging. Methods Firstly, using the improved complex nonlinear diffusion preprocessing filtering reduced the noise of the sequence images, then using the image mass center aligned the distortion data in the Y direction, finally, using the method of deviation average corrected motion artifacts along the Y direction. Results After correction, the motion artifacts in the longitudinal 2D images and the 3D image disappeared, the surface of the 2D and the 3D image became more smooth, the structure between layers of the images got clear and distinct, retinal en-face single slice image was sharp ,and the fundus tissue structure could be observed. Conclusion The algorithm of correction makes the physical structure of the retinal display truly after motion artifacts correction.

optical coherence tomography, motion artifacts, correction

TP391.41

A

1671-7104(2015)01-0005-04

10.3969/j.issn.1671-7104.2015.01.002

2014-09-12

天津科技创新体系及平台建设计划项目（13ZXGCCX06300）；中国医学科学院生物医学工程研究所院所基金(1417)

张明蓉，E-mail: zhmrong@yeah.net

李跃杰，E-mail: liyj_1@sina.com