时间序列三维视网膜SD-OCT图像的运动伪差校正方法与流程

本发明涉及出现处理技术领域，尤其涉及一种sd-oct图像的运动伪差校正方法。

背景技术：

视网膜是位于眼球后部的感光组织，是人类视觉系统的重要组成部分。oct是一种光学层析成像技术，其利用近红外弱光相干光照射样品组织，依据光的相干性产生干涉，对组织表层成像。具体，当oct(opticalcoherencetomography，光学相干断层扫描)系统对样品组织探测时，将样品组织等效为一系列的反射面，探测某深度上组织的反射率(即后向散射率)，得到其a扫描成像。将光束对组织样品进行横向移动，得到多个a扫描图像后进行拼接得到二维断层图像，即b扫描图像。如图1为视网膜三维sd-oct图像示意图，多个连续的a扫描(图示中a-scan)图像组成b扫描(图示中b-scan)图像，多幅连续的二维b扫描图像组成三维sd-oct成像。随着oct技术的迅速发展，sd-oct(spectral-domainopticalcoherencetomography，频域光学相干断层扫描)对于理解视网膜的结构和功能，解释视网膜疾病的发病机理，确定新型治疗方案，监测疾病治疗效果等方面起着越来越重要的作用，成为了眼科检查中的重要辅助手段，能快速、高分辨率的显示视网膜内部分层的三维图像。

sd-oct图像拍摄的过程一般需要持续几秒钟，在这段时间内，人的呼吸、心跳等不可避免的因素均会引起眼球的抖动，从而造成视网膜sd-oct图像的运动伪差。视网膜sd-oct图像的运动伪差主要有三种：轴向运动伪差，横断面运动伪差以及图像间运动伪差。如图2为视网膜三维sd-oct图像运动伪差示意图，其中，图2(a)为时间点1采集的sd-oct图像，图2(a1)为该sd-oct图像中y方向的视图(y-z视图)，图2(a2)为该sd-oct图像中x方向的视图(x-z视图)，图2(a3)为该sd-oct图像垂直方向的投影图像(x-y视图)；图2(b)为时间点2采集的sd-oct图像，图2(b1)为该sd-oct图像中y方向的视图(y-z视图)，图2(b2)为该sd-oct图像中x方向的视图(x-z视图)，图2(b3)为该sd-oct图像垂直方向的投影图像(x-y视图)。具体，由于眼睛和头部在垂直方向的运动，在sd-oct图像的慢扫描方向，即y方向(图2(a)中y-z视图间的伪差，图2(b)中y-z视图间的伪差)，相邻的两个b扫描图像之间往往存在很大的失真，主要表现为相邻的两个b扫描图像的跳变，称为轴向运动伪差。眼睛和头部在水平方向的运动也会引起失真，主要表现为sd-oct图像投影图中血管的错位，如图2中圈q所在的位置，称为横断面运动伪差。横断面运动伪差又可以分为两种，第一种为面内运动伪差，主要表现为快扫描方向，即x方向的图像扫描跳跃；第二种为面外运动伪差，也是最难校正的一种运动伪差，主要表现为同一位置重复扫描或者漏扫描。轴向运动伪差和横断面运动伪差都是sd-oct图像内部的运动伪差，但是，在进行眼科检查时，往往需要对一个时间段内的视网膜情况进行观察比较，因此，对于时间序列三维视网膜sd-oct图像来说，除了图像内部的运动伪差外，不同时间点的sd-oct图像之间由于拍摄的角度和位置的不同还存在图像间的运动伪差(图2(a)和图2(b)间的伪差)。三维视网膜sd-oct图像的运动伪差给后期图像的处理以及医生的准确诊断等带来了很大困难，因此对运动伪差进行校正是非常有必要的。

校正视网膜sd-oct图像运动伪差的方法有很多，主要可以分为两大类：基于硬件的方法和基于软件的方法。其中，基于硬件的方法通过在oct扫描仪上增加相关硬件来搜集眼动的信息，利用搜集的眼动信息对sd-oct图像运动伪差进行实时或后期校正处理；基于软件的方法通过将采集的图像与参照图像进行对比或者直接对眼睛的运动规律进行科学假设来获得眼动信息。在这两种方法中，由于基于硬件的方法可以直接获得真实的眼动信息，因此消除运动伪差的效果比基于软件的方法好，但是缺点在于需要增加额外的硬件，会浪费大量的资源。此外，很多以前拍摄sd-oct图像都没有直接的眼动信息，因此只能用通过软件的方法来解决运动伪差的问题。

在基于软件的方法中，常用的方法包括基于正交扫描图像的方法，该方法通过参照两幅正交的扫描图像来获得眼动信息实现运动伪差的校正。但是，在临床上，通常并不会特意采集正交扫描图，因此，该方法并不适用于临床上采集的视网膜sd-oct图像的运动伪差校正。

技术实现要素：

针对上述现有技术的不足，本发明提供了一种时间序列三维视网膜sd-oct图像的运动伪差校正方法，有效解决现有技术中不能有效对sd-oct图像的运动伪差进行校正的技术问题。

为了实现上述目的，本发明通过以下技术方案实现：

一种时间序列三维视网膜sd-oct图像的运动伪差校正方法，包括：

s11获取基于时间序列的多张三维视网膜sd-oct图像，并从中选定一sd-oct图像作为其他待校正sd-oct图像的基准图像；

s12采用多尺度三维图搜索方法依次对获取的各sd-oct图像中的视网膜内部进行分层，得到包括视网膜色素上皮层在内的多个分层；

s13提取视网膜色素上皮层的灰度数据，并根据该灰度数据得到包含血管投影信息的垂直投影图像；

s14采用基于海森矩阵的多尺度血管增强滤波器对垂直投影图像中血管的管状结构进行检测；

s15采用形态学处理方法对检测到的血管管状结构进行腐蚀，提取得到血管的骨架结构；

s16采用surf算法提取血管骨架结构的特征点；

s17基于提取的特征点采用刚性配准模型得到各待校正sd-oct图像与基准图像之间的变换参数，依次完成时间序列sd-oct图像中待校正sd-oct图像与基准图像间运动伪差的校正。

进一步优选地，在步骤s13中，提取视网膜色素上皮层的灰度数据之后，计算视网膜色素上皮层上表面和下表面间垂直方向上的灰度平均值，得到包含血管投影信息的垂直投影图像。

进一步优选地，在步骤s16中，采用surf算法提取血管骨架结构的特征点之后，还包括使用随机抽样一致性方法对提取的特征点进行优化选择得到匹配特征点的步骤；

在步骤s17中，基于提取的匹配特征点采用刚性配准模型得到各待校正sd-oct图像与基准图像之间的变换参数，依次完成时间序列sd-oct图像中待校正sd-oct图像与基准图像间运动伪差的校正。

进一步优选地，在步骤s13之后，还包括：

s18采用直方图均衡化方法对垂直投影图像进行图像增强处理；

s19采用维纳滤波器对图像增强后的垂直投影图像进行去噪。

进一步优选地，在步骤s12之后，还包括：

s21根据视网膜色素上皮层下表面中的最低处依次对各分层的表面进行平坦化，完成sd-oct图像的轴向运动伪差校正。

进一步优选地，完成了sd-oct图像间运动伪差的校正和轴向运动伪差校正之后，还包括：

s31依次将时间序列的sd-oct图像进行投影得到相应的二维投影图像；

s32计算待校正二维投影图像偏离平均误差的方差σi,j：

其中，i表示待校正二维投影图像采样的时间点，n表示待校正二维投影图像采样时间点的数量，0≤n≤i-1；sadi-n,j表示时间点i采样的二维投影图像与时间点n采样的二维投影图像第j行中各像素点绝对差值的总和；

s33根据计算得到的方差进一步计算方差σ'i,j：

其中，σi-1,j表示时间点i-1采集的待校正二维投影图像偏离平均误差的方差；

s34根据计算方差σ'i,j及预先设定的灵敏度参数ω进一步计算自适应运动伪差阈值γ：

s35判断方差σi,j是否大于自适应运动伪差阈值γ，若是，跳转至步骤s36；否则，判定时间点i采样的二维投影图像不存在横断面运动伪差；

s36判定时间点i采样的二维投影图像存在横断面运动伪差，并将时间点i采样的二维投影图像与时间点i-1采样的二维投影图像进行配准拼接，完成其横断面运动伪差的校正。

本发明提供的时间序列三维视网膜sd-oct图像的运动伪差校正方法，无需增加额外硬件设备、也无需额外采集正交扫描图像即可完成sd-oct图像的运动伪差校正，包括轴向运动伪差校正、横断面运动伪差校正及图像间运动伪差校正，方法简单，大量节约资源的同时得到精确的校正结果，对于患者的病情发展评估及医生的术前以及术后的治疗方案的确定起到了重要的辅助作用。

附图说明

结合附图，并通过参考下面的详细描述，将会更容易地对本发明有更完整的理解并且更容易地理解其伴随的优点和特征，其中：

图1为视网膜三维sd-oct图像示意图；

图2为视网膜三维d-oct图像运动伪差示意图；

图3为本发明中时间序列三维视网膜sd-oct图像的运动伪差校正方法一种实施方式流程示意图；

图4为本发明一实例中视网膜内部分层示意图；

图5为本发明中非刚性配准示意图；

图6为本发明一实例中图像间运动伪差校正过程示意图；

图7为本发明中时间序列三维视网膜sd-oct图像的运动伪差校正方法另一实施方式流程示意图；

图8为本发明一实例中轴向运动伪差校正过程示意图；

图9为本发明一实例中三维视网膜sd-oct图像运动伪差校正效果三维渲染图。

具体实施方式

为使本发明的内容更加清楚易懂，以下结合说明书附图，对本发明的内容作进一步说明。当然本发明并不局限于该具体实施例，本领域内的技术人员所熟知的一般替换也涵盖在本发明的保护范围内。

如图3所示为本发明提供的时间序列三维视网膜sd-oct图像的运动伪差校正方法一种实施方式流程示意图，从图中可以看出，在该运动伪差校正方法中包括：

s11获取基于时间序列的多张三维视网膜sd-oct图像，并从中选定一sd-oct图像作为其他待校正sd-oct图像的基准图像；

s12采用多尺度三维图搜索方法依次对获取的各sd-oct图像中的视网膜内部进行分层，得到包括视网膜色素上皮层在内的多个分层；

s13提取视网膜色素上皮层的灰度数据，并根据该灰度数据得到包含血管投影信息的垂直投影图像；

s14采用基于海森矩阵的多尺度血管增强滤波器对垂直投影图像中血管的管状结构进行检测；

s15采用形态学处理方法对检测到的血管管状结构进行腐蚀，提取得到血管的骨架结构；

s16采用surf算法提取血管骨架结构的特征点；

s17基于提取的特征点采用刚性配准模型得到各待校正sd-oct图像与基准图像之间的变换参数，依次完成时间序列sd-oct图像中待校正sd-oct图像与基准图像间运动伪差的校正。

在本实施方式中，基于时间序列的多张三维视网膜sd-oct图像具体由时间序列中不同时间点(时间点1、时间点2、时间点3…)采样的三维视网膜sd-oct图像组成。在对图像间运动伪差进行校正的过程中，首先选定一sd-oct图像作为基准图像，一般来说，选定该时间序列中第一个时间点(对应上述时间点1)采样的sd-oct图像作为基准图像，即时间序列下的第一张sd-oct图像，当然，在实际应用中也可以根据实际情况进行基准图像的选定，这里不做具体限定。在选定了基准图像之后，该时间序列下的其他sd-oct图像以基准图像为基准进行图像间运动伪差的校正。

视网膜内部分层对于分析视网膜病变，如眼外伤的严重程度、黄斑水肿的形成等具有重要作用，故在选定了基准图像之后，采用多尺度三维图搜索方法依次对获取的各sd-oct图像中的视网膜内部进行分层，得到包括视网膜色素上皮层在内的多个分层。多尺度三维图搜索方法基于图论的方法，采用从粗糙到精细的不同分辨率，对视网膜内部各表面进行检视；同时基于边界的代价函数，当代价函数最小时，即找到各个表面。

在一实例中，如图4所示，将视网膜内部自动分割成6层(对应7个表面)，包括：神经纤维层rnfl、神经节细胞层与内丛状层gcl+ipl、内核层inl、外丛状层opl、外核层与内节层onl+isl、视网膜色素上皮层rep+，由分层从上到下的顺序，神经纤维层rnfl的上表面为第1表面，神经节细胞层与内丛状层gcl+ipl的上表面为第2表面(同神经纤维层rnfl的下表面)，以此类推，视网膜色素上皮层rep+的上表面为第6表面、下表面为第7表面。在其他实例中，还可以采用其他方法对视网膜内部进行分层得到相应的层数，这里不做具体限定，只要分层的结果中包括视网膜色素上皮层rep+即包含在本实施方式的内容中。

完成了视网膜内部的分层之后，随即提取视网膜色素上皮层的灰度数据，并计算视网膜色素上皮层上表面和下表面间垂直方向上的灰度平均值，得到包含血管投影信息的垂直投影图像。在如图4中的实例中，即提取第6表面和第7表面之间的灰度值数据，并计算第6表面和第7表面间垂直方向上的灰度平均值得到包含血管投影信息的垂直投影图像。由于sd-oct图像散班噪声的影响，垂直投影图也含有噪声，所以为了提升垂直投影图像的质量，需要对其进行图像增强处理及去噪处理，具体，采用直方图均衡化方法对垂直投影图像进行图像增强处理之后，采用维纳滤波器对其进行去噪，得到去噪后的垂直投影图像。

血管的准确检测和提取是sd-oct图像配准的一个重要步骤，在对垂直投影图像进行图像增强和去噪之后，本实施方式中首先采用了基于海森(hessian)矩阵的多尺度血管增强滤波器对血管的管状结构进行检测；随后采用形态学处理方法对检测到的管状血管进行腐蚀，从而提取到血管的骨架结构；接着对提取到的血管骨架结构采用surf特征(加速稳健特征)算法提取特征点；之后采用随机抽样一致性方法(一种迭代算法，根据一组包含异常数据的样本数据集，计算出数据的数学模型参数，从而得到有效的样本数据)对特征点进行优化选择，得到能够进一步加快运算速度和减少错误的匹配特征点。最后，基于提取的匹配特征点采用刚性配准模型(非刚性配准如图5所示)得到各待校正sd-oct图像与基准图像之间的变换参数，如旋转参数、平移参数等，依次完成时间序列sd-oct图像中待校正sd-oct图像与基准图像间运动伪差的校正。

如图6所示为一实例中图像间运动伪差校正过程示意图，其中，图6(a)为时间点1拍摄的sd-oct图像的垂直投影图像中提取的血管骨架结构(基准图像)，图6(b)为时间点2拍摄的sd-oct图像的垂直投影图像中提取的血管骨架结构(待校正sd-oct图像)，图6(c)为基准图像和待校正sd-oct图像中血管骨架结构的匹配特征点匹配结果，图6(d)为采用刚性配准模型配准前基准图像和待校正sd-oct图像的血管骨架结构重叠示意图，图6(e)为采用刚性配准模型配准后基准图像和待校正sd-oct图像的血管骨架结构重叠示意图。从中可以看出，配准之后，大大降低了基准图像和待校正sd-oct图像间的运动伪差。

对上述实施方式进行改进得到本实施方式，在本实施方式中，包括上述对图像间运动伪差进行校正的步骤之外，还包括轴向运动伪差校正的步骤，具体，如图7所示，在步骤s12，采用多尺度三维图搜索方法依次对获取的各sd-oct图像中的视网膜内部进行分层，得到包括视网膜色素上皮层在内的多个分层之后，还包括：s21根据视网膜色素上皮层下表面中的最低处(鲁棒性相对较高)依次对各分层的表面进行平坦化，完成sd-oct图像的轴向运动伪差校正。参照图4所示的分层结果，这里即参考第7表面最低处位置对其他分层进行平坦化，平坦化后的sd-oct图像如附图8所示，其中，图8(a)为校正前sd-oct图像示意图，图8(b)为校正后sd-oct图像示意图。

对上述实施方式进行改进得到本实施方式，在本实施方式中，完成上述图像间运动伪差的校正和轴向运动伪差校正之后，还包括横断面运动伪差校正步骤，完成sd-oct图像横断面运动伪差的校正，尤其是当某个时间点的sd-oct图像出现漏扫描时，借助其他时间点的采样sd-oct图像对漏扫描的信息进行补偿。

在横断面运动伪差校正的过程中，首先，依次将完成了图像间运动伪差和轴向运动伪差校正的sd-oct图像进行投影得到相应的二维投影图像(对应b扫描图像)；之后，参照式(1)计算待校正二维投影图像偏离平均误差的方差σi,j：

其中，i表示待校正二维投影图像采样的时间点，n表示待校正二维投影图像采样时间点的数量，0≤n≤n-1；sadi-n,j表示时间点i采样的二维投影图像与时间点n采样的二维投影图像第j行中各像素点绝对差值的总和。

之后，根据式(2)计算得到的方差进一步计算方差σ'i,j：

其中，σi-1,j表示时间点i-1采集的待校正二维投影图像偏离平均误差的方差。

接着，根据计算方差σ'i,j及预先设定的灵敏度参数ω进一步计算自适应运动伪差阈值γ，如式(3)：

其中，ω为控制横断面运动伪差识别灵敏度，ω越小，自适应运动伪差阈值对横断面运动伪差越敏感，在实际应用中，ω的值根据实际情况进行选定，如将其选定为0.5、1等。

最后，判断方差σi,j是否大于自适应运动伪差阈值γ，若是，判定时间点i采样的二维投影图像存在横断面运动伪差，记录其位置，并将时间点i采样的二维投影图像与时间点i-1采样的二维投影图像进行配准拼接，完成其横断面运动伪差的校正。否则，判定时间点i采样的二维投影图像不存在横断面运动伪差。二维投影图像进行配准拼接的过程具体包括：在运动伪差处将相对应的二维投影图像与前一个采样点的二维投影图像进行基于图像灰度值的配准，找到最佳平移参数并利用该平移参数对横断面运动伪差进行校正。之后，将校正后的图像再一次利用公式(1)～(3)进行横断面运动伪差检测，若无伪差，说明该处伪差为横断面内运动伪差，且已校正完毕。若还有伪差存在，则判断为横断面运动外伪差，并借助时间序列的sd-oct图像对漏扫描图像进行填补。具体，在横断面运动伪差处将二维投影图像分裂为子图像，以最大的一块子图像为基础图像，其他子图像包括前一个时间点的图像均根据surf特征法(提取surf特征点)进行图像配准，找到匹配位置，然后进行图像融合，以克服重复扫描或者获得漏扫描的信息，达到修正横断面运动外伪差的目的。

图9为一实例中三维视网膜sd-oct图像运动伪差校正效果三维渲染图，其中，图9(a1)和图9(b1)为校正前时间点1和时间点2的三维视网膜sd-oct图像，图9(a2)和图9(b2)为校正后时间点1和时间点2的三维视网膜sd-oct图像。