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

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

（1.天津大学精密仪器与光电子工程学院，天津 300072；2.上海美沃精密仪器股份有限公司，上海 200237）

1 引言

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

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

2 OCT 图像错位分析

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

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

图1 OCT 序列图像错位原因。（a）眼睛轴向运动；（b）眼睛垂轴运动；（c）眼睛旋转运动；（d）多角度光线入射Fig.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

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

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

3 算法设计

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

图2 配准算法流程图Fig.2 Flow chart of the registration algorithm

3.1 预处理

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

图3 预处理过程。（a）原始视网膜B-Scan 图像；（b）列平均灰度分布图；（c）ROI 图像Fig.3 Pretreatment process.(a)Original B-Scan image of retinal image;(b)column average gray distribution map;(c)ROI image

3.2 平移配准

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

POC 算法[22−24]利用基于傅立叶变换的平移不变性将图像在空域中的平移量转化成对POC函数峰值位置的求取。假设目标图像f2(x,y)是由参考图像f1(x,y)平移得到，对二者进行离散傅立叶变换得到F1(u,v)和F2(u,v)，计算二者互功率谱的相位G(u,v)，再进行反傅立叶变换可得到POC 函数：

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

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

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

图4 分段拟合处理过程。（a）曲线分割结果（黑色虚线为分割线）；（b）原始曲线段（最佳连接点A、黑色直线为包络线）；（c）处理后的曲线段Fig.4 Fitting process of curve segments.(a)Results of curve segmentation(the dotted black line is the segmentation line);(b)original curve segment results(optimal connection point A,the straight black line is the envelope);(c)curve segment after processing

3.3 旋转配准

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

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

其中n为总点数，根据统计学规律可知，当满足|ym+1−ym|<0.95 md条件时剔除错误点，并采用双线性插值法得到完整、准确的特征点集合。图5（a）（彩图见期刊电子版）是上边界提取过程的示意图，A点为下边缘区域的最大灰度像素点，红色单箭头线为行局部搜索范围，红色曲线为拟合后的视网膜上边界。

图5 刚性配准结果。（a）上边界提取过程；（b）平移配准后的叠加图；（c）旋转配准后的叠加图（绿色为参考图像，紫色为目标图像）Fig.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)

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

3.4 畸变校正

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

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

其中，Io(x,y+∆y)是目标图像在坐标(x,y+∆y)处的灰度值，Ir(x,y)是参考图像在坐标(x,y)处的灰度值，D为图像对上边界之间的相对距离。能量函数反映得是一对一轴向扫描的灰度分布差异，在一定的轴向平移范围(−2D≤∆y≤2D)内寻找E(∆y)的最小值，得到高精度的轴向匹配结果（图6，彩图见期刊电子版）。非刚性配准后的叠加图及局部放大图如图6（b）所示，可见，相较于刚性配准后的叠加图及局部放大图（6（a）），非刚性配准校正了图像的边缘畸变，避免了重影区域的产生，提高了图像的配准精度。

图6 非刚性配准结果。（a）刚性配准后的叠加图及局部放大图；（b）非刚性配准后的叠加图及局部放大图（绿色为参考图像，紫色为目标图像）Fig.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)

4 测量实验与结果

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

为了评价本文方法的性能，采用信噪比（Signal-to-Noise Ratio,SNR）和对比度（Contrast-to-Noise Ratio,CNR）的客观标准来定量分析图像质量，相应的计算公式[27]分别为：

式中，µb和 σb分别为背景区域的平均值和标准差，µe和 σe分别为第k个视网膜有效区域的平均值和标准差，K为有效区域的总个数。以图8（a）（彩图见期刊电子版）所示的原始单帧图像为例，选定一个较大的黄色矩形框为背景区域和7 个红色矩形框为视网膜有效区域（K=10）。利用本文方法、文献[18]方法和文献[21]方法对多帧连续图像进行配准，计算不同帧数叠加图在规定区域下的SNR 和CNR 值，相应的变化趋势分别如图8（b）（彩图见期刊电子版）和图8（c）（彩图见期刊电子版）所示。

图7 本算法实验结果图。（a）原始单帧图；（b）配准后的叠加图Fig.7 The experimental results obtained by the algorithm proposed in this paper.(a)Original single frame image;(b)averaging image after registration

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

图8 （a）区域选择及SNR（b）和CNR（c）随叠加帧数递增的曲线变化图Fig.8 (a)Area selection,SNR(b)and CNR(c)curve change with an increasing number of multiple frames

表1 本文配准算法较对比算法的提升效果Tab.1 Improvement performance of proposed algorithm comparing with other algorithms

5 结论

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