CN106447708A - 一种oct眼底图像数据配准方法 - Google Patents

Abstract

本发明选取OCT眼底图像，使用Canny边缘检测法提取OCT眼底图像的视网膜边缘，并将这些边缘以点云数据的格式收集起来；其次，采用空间网格划分的方法提取点云数据的特征点；再次，应用奇异值分解算法计算出待配准点云间的变换矩阵，消除明显的位置误差；最后，使用改进的迭代最近点算法进行精确配准，把得到的旋转矩阵和平移矩阵应用于原始的OCT眼底图像上，得到最终的结果。当处理具有较大数据量的密集型点云时，本发明算法在时间复杂度和配准精度方面有明显的优势。在大多数情况下，本发明算法将传统的迭代最近点算法效率提高了70％，本发明不仅对OCT眼底图像配准与拼接有效，并且创建具有较大视野眼底视网膜图像精确性。

Description

一种OCT眼底图像数据配准方法

技术领域

本申请涉及图像配准技术，具体的，涉及对于采用OCT成像的三维眼底图像数据进行配准的方法。

背景技术

图像配准技术通常是指对于一幅图像，寻找某一系列的空间变换，使其与另一幅或多幅图像对应的特征信息有相同的空间位置。对于OCT三维成像的眼底数据而言，配准的结果应使不同眼底图像中视网膜各个层相互对齐而不会出现明显的断层。

目前图像配准技术主要分为两类：基于特征的图像配准和基于互信息的图像配准。基于特征的图像配准方法首先提取图像的特征信息，然后以这些特征为模型进行配准。特征提取的结果是一组包含图像特征的数据以及对这些数据的描述，每个数据由一组属性表示，对属性的进一步描述包括边缘的定向和弧度的大小等。这些特征构成了图像的局部特征，而局部特征之间存在着相互关系，如几何关系、辐射度量关系、拓扑关系等。可以用这些局部特征之间的关系描述全局特征。通常基于局部特征配准大多都是基于点、线或边缘的，而全局特征的配准则是利用局部特征之间的关系进行配准的方法。基于互信息的配准方法通常是指最大互信息的图像配准方法。基于互信息的图像配准是用两幅图像的联合概率分布与完全独立时的概率分布的广义距离来估计互信息，并作为多模态医学图像配准的测度。当两幅基于共同场景的图像达到最佳配准时，它们的对应像素的灰度互信息应为最大。由于基于互信息的配准对噪声比较敏感，一般情况下，需要通过滤波和分割等方法对图像进行预处理，接着进行采样、变换、插值、优化等过程达到配准的目的。

二维医学图像配准的文献有很多，Kratika Sharma和Ajay Goyal(Sharma K,Goyal A.Classification based survey of image registrationmethods.International Conference on Computing,Communications&NetworkingTechnologies,2013,1-7)将图像配准方法步骤分为四步，分别为空间关系方法，放缩方法、金字塔方法和恒定描述符小波方法。Mei-sen Pan等学者(Pan M,Jiang J,Rong Q,et al.Amodified medical image registration.Multimedia Tools&Applications,2013,70:1585-1615)提出了一种基于边缘检测的B样条梯度算子方法。这种方法实现简单，并具有较低的计算负载和良好的配准精度。而对于两个完全重叠且存在噪声的图像，LucianCiobanu和Luís(Ciobanu L,L.Iterative filtering of SIFTkeypoint matches for multi-view registration in Distributed VideoCoding.Multimedia Tools&Applications,2011,55:557-578)采用迭代的策略重建点对间的关系并以此来消除噪声点。这种策略使得噪声点的总体数量显著减少，同时又保持了原有的高速精确匹配。为了使图像配准算法更加准确和健壮，有一种新颖的想法是建模整张图像分布模型。Shihui Ying等学者(Ying S,Wu G,Wang Q,et al.GroupwiseRegistration via Graph Shrinkage on the Image Manifold.IEEE Conference onComputer Vision&Pattern Recognition,2013,2323-2330)最先引入这个概念，并将图像特征点配准的步骤可以归纳为动态收缩问题。除上述方法外，其他一些基于特征的配准算法也能得到较好的结果。

然而，光学相干性断层扫描眼底数据是由多个二维光学相干性断层扫描图像叠加组成的三维体数据。以上方法在处理这样的数据时将面临运行内存和计算时间的限制。因此，为了配准三维眼底数据，一个得到具有较大视野光学相干性断层扫描体数据的方案是采用蒙太奇方法(Li Y,Gregori G,Lam B L,et al.Automatic montage of SD-OCT datasets.Optics express,2011,19:26239-26248)。此方法使用以血管脊为兴趣特征，利用重采样、插值和互相关的方法来拼凑完整的光学相干性断层扫描体数据。这种蒙太奇方法可以将分散，部分重叠的光学相干性断层扫描图像拼接成一个具有较大视野的三维图像。然而，当眼底图像中作为兴趣特征的血管脊模糊时，这种方法将无法完成配准。除此之外，也有学者利用现有的工具和平台来生成具有较大视野体数据。Meng Lu(MengL.Acceleration method of 3D medical images registration based on computeunified device architecture.Bio-medical materials and engineering,2014,24:1109-1116)基于由NVIDIA提供的计算设备架构出了一个加速配准的运算方法，该算法可以提高三维医学图像配准的性能，加快计算速度，适合处理大规模数据。此外，StephanPreibisch等学者(Preibisch S,Saalfeld S,Tomancak P.Globally optimal stitchingof tiled 3D microscopic image acquisitions.Bioinformatics,2009,25:1463-1465)还在imageJ平台上实现了一个拼接插件，它能够在不需要先验知识的情况下将分散的小块体数据重建成一个整体。除上述拼接插件外，其他类型的拼接工具也陆续开始应用。然而，由于光学相干性断层扫描设备自身的瓶颈以及扫描过程中人眼的不自觉的运动，所得到的光学相干性断层扫描眼底数据会掺杂一些微小的非刚性变换(陈国琳.非刚性医学图像配准方法的研究与实现[D].南京理工大学硕士论文,2009)。因此上述方法在处理此类较为特殊的临床眼科光学相干性断层扫描图像时具有一定的局限性。

综上，对于OCT三维成像的眼底数据，如何能够快速、精确的创建具有较大视野的眼底图像的算法，以此为临床医师对眼科疾病的诊断和治疗提供帮助，是现有技术亟需解决的技术问题。

发明内容

本发明的目的在于提出一种OCT眼底图像数据配准方法，以快速、精确创建具有较大视野范围的眼底视网膜图像数据。

为达此目的，本发明采用以下技术方案：

一种OCT眼底图像数据配准方法，包括如下步骤：

步骤一、采用Canny边缘检测方法对图像进行去噪、边缘检测：

对原始图像进行灰度化，对图像进行滤波，然后计算图像的梯度幅值，对梯度幅值进行非极大值抑制，然后用双阈值算法检测和连接边缘；

步骤二、图像点云数据的提取及可视化：

将经过步骤一处理的图像按照原有的顺序叠加成一个三维体数据，并把每一个边缘点按照其所处的空间位置抽取成三维点云中的一个点，点云点的三维坐标对应叠加得到的三维体数据的空间坐标；

步骤三、点云数据边界特征点的提取：

将点云中的点按照其空间坐标分配到不同的空间网格中，接着找出所有属于点云边界的空间网格，最后将边界网格中的点云点提取出来作为点云数据的边界特征点；

步骤四、采用奇异值分解算法完成点云数据初始配准：

令P代表原始集合，Q代表对照集合，参见公式3，定义目标矩阵如下：

C_P和C_Q分别是原始集合P以及对照集合Q的质心，M代表点云数据集中点云的数量，P_i和Q_i分别代表原始集合P与对照集Q合中第i个点，将目标矩阵E采用奇异值分解(SingularValue Decomposition,SVD)，则E＝UDV^T。其中U的列是EE^T矩阵的特征向量，V的列是E^TE矩阵的特征向量。V^T是V的转置矩阵且D为对角矩阵，D＝diag(d_i)，其中，d_i为E的奇异值，令

则旋转矩阵R＝UBV^T，平移矩阵T＝C_Q-RC_P，将所求得的矩阵R和T作用于原始集合P以消除初始条件下原始集合P与对照集合Q可能存在的较大位移误差；

步骤五、采用改进迭代最近点算法完成点云数据精确配准：

步骤5.1初始化。对于给定的两个点云集合P和Q，指定一个收敛阈值τ，

步骤5.2采用公式5算出点云数据中每个点对应的权值，比较阈值ε并将权值较低的点排除，

Q_B为集合Q中点P_A的对应点，Dis(P_A,Q_B)代表P_A和Q_B间的欧氏距离，Dis_MAX代表点对间欧氏距离的最大值，欧式距离采用公式6计算得到，

步骤5.3迭代以下步骤，直到公式7的最小均方根误差收敛于所给定的阈值τ：

步骤5.3.1根据公式6计算出集合P和Q中点云之间的欧式距离，

其中，ω_x，ω_y，ω_z分别代表M-估计在各个坐标方向的权重因素，(x_A,y_A,z_A)，(x_B,y_B,z_B)分别为集合P中点A和Q中点B的空间坐标，

步骤5.3.2对于排除了权值较低的点的集合P，在集合Q中寻找欧氏距离最近的点作为对应点并存储在最近点集中，

步骤5.3.3利用公式7采用最小二乘法计算计算集合P和最近点集间的旋转矩阵R和平移矩阵T，

步骤5.3.4将旋转矩阵R和平移矩阵T应用于集合P，得到新的集合，利用公式7计算最小均方根误差是否收敛于所给定的阈值τ，如果是则结束运算，否则利用步骤5.3进行迭代计算。

进一步的，在步骤一中，原始图像的灰度化公式为Gray＝0.299R+0.587G+0.114B，图像滤波采用高斯滤波，采用一阶有限差分来计算图像像素矩阵关于横向和纵向的偏导数。

进一步的，在步骤一中，采用双阈值来检测并连接边缘，高阈值的作用在于抑制边缘，所有高于此阈值的梯度幅值才能被视为边缘，低阈值的作用在于连接边缘，将所有梯度幅值高于低阈值的点视为边缘并连接起来成为最终边缘检测的结果。

进一步的，在步骤三中，采用有向包围盒作为点云数据的最小包围盒，以将点云数据点划分到不同的空间网格中，

按照相同的体积将点云数据分解成若干个空间网格，每个空间网格的大小被定义为：

公式中L代表点云数据中点云点的数量，V代表最小包围盒的体积。V/L为点云密度的倒数，代表点云数据中每一个点云点所占空间的平均大小。令空间网格的初始大小S_grid为的点云密度倒数的K倍，并将点云数据最小包围盒按照S_grid大小分成若干个空间网格。

进一步的，使用边界种子网格算法找到所有边界空间网格，并提取这个边界空间网格中包含的点云点作为点云数据的边界特征点，将空间网格分成两类，空格和实格，不包含任何点云点的空间网格为空格，其他的空间网格为实格。以空间坐标(x,y,z)表示某一网格空间位置，用函数f表示此空间网格的类型，如果某个空间网格是实格，f(x,y,z)＝1，否则f(x,y,z)＝0，使用公式2来判定某个空间网格是否为的边界空间网格：

当U(x,y,z)≤1时，代表此空间网格的六邻域中上下、左右、前后的网格中至多有一对都是实格，那么此空间网格是一个边界网格。

进一步的，在步骤三中，8≤K≤24。

进一步的，在步骤五中，收敛阈值τ＝0.2，比较阈值ε为，0.2≤ε≤0.4。

进一步的，在步骤五中，M-估计在各个坐标方向的权重公式为

v为各个坐标方向上标准化的残差值，c为常数。

进一步的，c＝1.345。

本发明选取OCT眼底图像，使用Canny边缘检测法提取OCT眼底图像的视网膜边缘，并将这些边缘以点云数据的格式收集起来；其次，采用空间网格划分的方法提取点云数据的特征点；再次，应用奇异值分解算法计算出待配准点云间的变换矩阵，消除明显的位置误差；最后，使用改进的迭代最近点算法进行精确配准，把得到的旋转矩阵和平移矩阵应用于原始的OCT眼底图像上，得到最终的结果。

当处理具有较大数据量的密集型点云时，本发明算法在时间复杂度和配准精度方面有明显的优势。在大多数情况下，本发明算法将传统的迭代最近点算法效率提高了70％。本发明不仅对OCT眼底图像配准与拼接有效，并且创建具有较大视野眼底视网膜图像精确性。

附图说明

图1是根据本发明的具体实施例的OCT眼底图像数据配准方法的流程图；

图2是根据本发明的具体实施例的采用边界网格种子算法的边界网格的状态图；

图3是根据本发明的具体实施例的点云数据配准过程对比图；

图4是根据本发明的具体实施例的点云数据重合区域部分结果放大对比图；

图5是根据本发明的具体实施例的变换矩阵作用于原始集合与对照集合的配准后位置对比图；

图6是根据本发明的具体实施例的对于“斯坦福兔子”点云数据的配准结果的对比图；

图7是根据本发明的改进的算法与传统迭代算法时间消耗的对比图。

具体实施方式

下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部结构。

本发明的原理为：采用Canny边缘检测方法对医学眼底图像进行去噪，并将得到的边缘作为图像的特征，将这些特征提取成点云数据，使用最小包围盒的方法提取出点云数据中的边界特征点，利用边界特征点采用奇异值分解完成初始配准，通过改进迭代最近点方法对点云数据进行精确配准。

参见图1，示出了根据本发明的OCT眼底图像数据配准方法的流程图，具体包括如下步骤：

步骤一、采用Canny边缘检测方法对图像进行去噪、边缘检测：

图像的边缘指图像从某一灰度值急剧变化到另一灰度的临界区域，它通常分隔开两快亮度变化显著的区域。图像的边缘携带了图像大部分信息，是图像特征提取的重要部分，该步骤一中包括如下子步骤：

(1)对原始图像进行灰度化：实际应用场景中，大部分图像为彩色格式。首先需要将彩色格式的图像灰度化。较为常见的灰度化方法为对图像RGB三个通道的采样值进行加权平均，即Gray＝(R+G+B)/3。对于医学眼科图像，考虑到人眼的生理特点，本发明中灰度化公式为Gray＝0.299R+0.587G+0.114B。

(2)对图像进行滤波，然后计算图像的梯度幅值：

图像滤波是图像处理中一个重要的过程，它在尽量保留图像原有特征的情况下去抑制和剔除噪声对图像处理的影响，从而提高图像处理的准确率与可靠性。常见的滤波方法主要包括均值滤波、中值滤波、高斯滤波和双边滤波等。本发明中Canny边缘检测算法采用高斯滤波。

图像的边缘一般处于灰度值变化较大的区域，可以使用一阶有限差分来计算图像像素矩阵关于横向和纵向的偏导数。对于偏导数越大的点，说明该点所处区域的灰度变化较大，有可能是图像的边缘点。

(3)对梯度幅值进行非极大值抑制，然后用双阈值算法检测和连接边缘：

由低阶导数的结论可知，梯度幅值为极大值的点不一定为区域灰度值变化最大的点，并不能说明此点一定为边缘点，所以本发明中Canny边缘检测算法进行了非极大值抑制。即搜寻梯度幅值极大值点的八邻域，将梯度幅值最大的点保留，其他的点略掉，确保任意小区域内最多不超过一个边缘点。

Canny边缘检测算法采用双阈值来检测并连接边缘。高阈值的作用在于抑制边缘，所有高于此阈值的梯度幅值才能被视为边缘。故高阈值越高，检测出来的边缘越少。低阈值的作用在于连接边缘，由于高阈值的抑制作用，导致许多检测出来的边缘点并不能连接成一条边。此时Canny边缘检测算法将边缘点中所有梯度幅值高于低阈值的点视为边缘并连接起来成为最终边缘检测的结果。

步骤二、图像点云数据的提取及可视化

点云数据由一组三维空间中的点构成，通常用来描述物体的三维空间结构。每一个点云点都至少包含一组表示自身空间位置的三维坐标，有些点云也具有高维坐标，这些坐标用来表示颜色或光照反射强度等信息。将原始眼底光学相干性断层扫描图像经过步骤一的处理后，得到了一系列包含视网膜各层边缘的图像。

将这些图像按照原有的顺序叠加成一个三维体数据，并把每一个边缘点按照其所处的空间位置抽取成三维点云中的一个点，点云点的三维坐标对应叠加得到的三维体数据的空间坐标。

步骤三、点云数据边界特征点的提取

在该步骤中，首先将点云中的点按照其空间坐标分配到不同的空间网格中，接着找出所有属于点云边界的空间网格，最后将边界网格中的点云点提取出来作为点云数据的边界特征点。

进一步的，要将点云数据点划分到不同的空间网格中，首先需要找到能够包含整个点云的最小包围盒，再将此最小包围盒等分成若干个空间网格，使得任意点云数据点都被包含在某个空间网格中。

本发明采用有向包围盒作为点云数据的最小包围盒，有向包围盒是一种方向任意的包围盒方法，可以更加紧凑的包含整个物体。得到点云数据的最小包围盒之后，本发明按照相同的体积将点云数据分解成若干个空间网格，每个空间网格的大小被定义为：

公式中L代表点云数据中点云点的数量，V代表最小包围盒的体积。V/L为点云密度的倒数，代表点云数据中每一个点云点所占空间的平均大小。令空间网格的初始大小S_grid为的点云密度倒数的K倍，并将点云数据最小包围盒按照S_grid大小分成若干个空间网格。公式1中K是一个可变参数，对于本发明所处理的密集型点云数据，实验证明当K被赋值在8～24时，空间网格的大小能够包含足够多的点云数据点。

进一步的，本发明使用边界种子网格算法找到所有边界空间网格，并提取这个边界空间网格中包含的点云点作为点云数据的边界特征点。为了提取边界空间网格，可以将空间网格分成两类，空格和实格，不包含任何点云点的空间网格为空格，其他的空间网格为实格。以空间坐标(x,y,z)表示某一网格空间位置，用函数f表示此空间网格的类型，如果某个空间网格是实格，f(x,y,z)＝1，否则f(x,y,z)＝0。空间中与此网格相邻的六邻域网格坐标分别为(x-1,y,z)，(x+1,y,z)，(x,y-1,z)，(x,y+1,z)，(x,y,z-1)，(x,y,z+1)。使用公式2来判定某个空间网格是否为的边界空间网格：

U(x,y,z)为三个乘积的和，令中间网格的坐标为(x,y,z)，那么(x-1,y,z)，(x+1,y,z)，(x,y-1,z)，(x,y+1,z)，(x,y,z-1)，(x,y,z+1)表示此网格的六邻域网格所在位置。f(x-1,y,z)·f(x+1,y,z)＝1当且仅当网格上下位置都是实格，f(x,y-1,z)·f(x,y+1,z)＝1当且仅当网格左右位置都是实格，f(x,y,z-1)·f(x,y,z+1)＝1当且仅当网格前后位置都是实格。当U(x,y,z)≤1时，代表此空间网格的六邻域中上下、左右、前后的网格中至多有一对都是实格，那么此空间网格是一个边界网格。

参见图2，是根据本发明的具体实施例的采用边界网格种子算法的边界网格的状态图，左图为平面边界，中间图为线边界，右图有点边界。

步骤四、采用奇异值分解算法完成点云数据初始配准

对于上述点云数据边界特征点，使用奇异值分解算法来计算待配准集合之间的旋转矩阵和平移矩阵，完成点云数据的初始配准。

令P代表原始集合，Q代表对照集合，参见公式3，定义目标矩阵如下：

C_P和C_Q分别是原始集合P以及对照集合Q的质心，M代表点云数据集中点云的数量，P_i和Q_i分别代表原始集合P与对照集Q合中第i个点，将目标矩阵E采用奇异值分解(SingularValue Decomposition,SVD)，则E＝UDV^T。其中U的列是EE^T矩阵的特征向量，V的列是E^TE矩阵的特征向量。V^T是V的转置矩阵且D为对角矩阵，D＝diag(d_i)，这里d_i为E的奇异值，即E'E矩阵的特征值的平方根，令

则旋转矩阵R＝UBV^T，平移矩阵T＝C_Q-RC_P，将所求得的矩阵R和T作用于原始集合P以消除初始条件下原始集合P与对照集合Q可能存在的较大位移误差，例如，P'＝RP+T，为精确配准步骤提供一个良好的状态。

步骤五、采用改进迭代最近点算法完成点云数据精确配准

对于传统的迭代最近点算法的主体思想是首先根据某种几何特征来寻找出P和Q的所有对应点，并将这些对应点作为配准的对象。接着构建表示配准程度的目标方程，通过寻找目标方程的最优解来计算出当前情况下的旋转矩阵和平移矩阵。最后不断地迭代此流程，直到目标方程满足某个设定的阈值。

传统的迭代最近点算法以点到点的距离作为几何特征来寻找对应点。对于集合P和Q中的两个点P_i(x_p,y_p,z_p)和Q_i(x_q,y_q,z_q)，它们之间的欧氏距离表示为

对于集合P中的任意一点，利用上述公式计算Q中每一点到此点的距离，并选取欧氏距离最小的点作为此点在集合Q中的对应点。接着寻找旋转矩阵R与平移矩阵T，将其作用于P_i，则所得到点的位置为RP_i+T，利用最小二乘法构造目标方程

其中，N代表点云中点的数量，E代表经过变换后集合P中每一个点到集合Q中对应点距离的平方和。

从上述目标方程中可以看出，当E最小时代表此次迭代对应点的相对位置最近。故使得E最小的旋转矩阵和平移矩阵便为此次迭代的最优解。对于旋转矩阵和平移矩阵的求解，采用平移和旋转分离的办法，即先对平移矩阵T进行初值估算,首先计算P和Q的重心分别为

和

则集合P和Q间的平移估算值为p_c-q_c，此时目标方程化简为

以此求得此次迭代过程中的旋转矩阵R并利用T＝Q-RP计算出平移矩阵T，重复迭代以上过程，直至目标方程最优解E收敛与给定的阈值。

对于具有较好初始条件的点云数据，迭代最近点算法能够获得较为精确的配准结果。然而，传统的迭代最近点算法也存在一些不足之处。

首先，迭代最近点算法假设两个点云集合为包含关系，即一个点云集合是另一个点云集合的子集，这一情况通常较难满足。

其次，在选取对应点对的过程中，对于点云数据集合中的任意一点，算法会计算此点与另一集合中所有点的欧氏距离。假设两个点云集合各有N_P和N_Q个点，那么此步骤的时间复杂度为O(N_P×N_Q)。这使得算法花费大量时间计算对应点对的欧氏距离，计算代价很大。此外，算法默认欧氏距离最近的点对为对应点，由于存在一定量的噪声点，故此步骤后产生一定的误差，使算法陷入局部最小值。

因此，本发明提出了改进迭代最近点算法：传统迭代最近点算法中，点云集合中的所有点被赋予了相同的权重，因此所有点都会参加到点与点间欧氏距离的计算过程中，这是算法的瓶颈所在。本发明算法对不同的点赋予不同的权重，点对间的欧氏距离越远，它们的权值就越小。假设P_A是集合P中的一点，则P_A的权重如公式5表示为

Q_B为集合Q中点P_A的对应点，Dis(P_A,Q_B)代表P_A和Q_B间的欧氏距离。Dis_MAX代表点对间欧氏距离的最大值,对于给定的阈值ε，所有权重低于ε的对应点对将被舍弃，不参与到迭代计算的过程中。阈值ε是一个可变参数，它用来权衡配准的时间复杂度和配准的精确度。若令ε减小，则意味着更多的对应点对会参与到迭代计算过程中，使配准结果更加精确，但与此同时也增加了计算次数，提高了算法的时间复杂度。优选地，ε为0.2～0.4。

由于点云数据可能会存在噪声，这些噪声会影响点云重心及迭代最近点目标方程的计算等过程，从而降低配准精度。本发明引入了M-估计以排除噪声点对实验结果的影响。M-估计稳健回归的基本思想是根据回归残差的大小确定各点的权重，以达到稳健的目的。为减少异常点的影响，可以对不同的点赋予不同的权重，即对残差小的点给予较大的权重，而对残差较大的点给予较小的权重。本发明中M-估计在各个坐标方向的权重公式为

某一点的权重因素被定义为

v为各个坐标方向上标准化的残差值，c为常数，一般取1.345。因此，当v属于区间(-c,c)时，M-估计退化为经典最小二乘估计。而当残差v大于c时，权重因素w_v随着残差的增大而减小。对应集合P中点A和Q中点B间欧氏距离的计算方法如公式6所示：

其中，ω_x，ω_y，ω_z分别代表M-估计在各个坐标方向的权重因素，(x_A,y_A,z_A)，(x_B,y_B,z_B)分别为点A和B的空间坐标，此时迭代方程如下：

综上，改进后的迭代最近点算法如下：

步骤5.1初始化。对于给定的两个点云集合P和Q，指定一个收敛阈值τ，优选地，τ为0.2；

步骤5.2采用公式5算出点云数据中每个点对应的权值，比较阈值ε并将权值较低的点排除，优选地，ε为0.2～0.4；

步骤5.3迭代以下步骤，直到公式7的最小均方根误差收敛于所给定的阈值τ：

步骤5.3.1根据公式6计算出集合P和Q中点云之间的欧式距离，

步骤5.3.2对于排除了权值较低的点的集合P，在集合Q中寻找欧氏距离最近的点作为对应点并存储在最近点集中，

步骤5.3.3利用公式7采用最小二乘法计算计算集合P和最近点集间的旋转矩阵R和平移矩阵T，

步骤5.3.4将旋转矩阵R和平移矩阵T应用于集合P，得到新的集合，利用公式7计算最小均方根误差是否收敛于所给定的阈值τ，如果是则结束运算，否则利用步骤5.3进行迭代计算。

因此，本发明选取OCT眼底图像，使用Canny边缘检测法提取OCT眼底图像的视网膜边缘，并将这些边缘以点云数据的格式收集起来；其次，采用空间网格划分的方法提取点云数据的特征点；再次，应用奇异值分解算法计算出待配准点云间的变换矩阵，消除明显的位置误差；最后，使用改进的迭代最近点算法进行精确配准，把得到的旋转矩阵和平移矩阵应用于原始的OCT眼底图像上，得到最终的结果。

实施例1：

下面结合实验，对本发明做进一步阐述，并将本发明和经典迭代最近点方法进行比较，验证其精确性和鲁棒性。

实验所使用的计算机配置为Intel Core E7500双核CPU，主频2.93GHz，内存为1GB×2 DDR2，操作***为Microsoft Windows Win7 Sp1旗舰版，算法的实现平台为MicrosoftVisual Studio 2010.

为了评估算法的配准性能，本发明算法采用两个眼底光学断层扫描图像作为实验数据，它们是人类视网膜结构相邻的部分。两个数据集合的重叠部分大约为75×500×375个体素。图3从四个不同的角度展示了点云数据的配准过程。图片的左方区域显示了两个点云集的初始位置(红色的点代表原始集合，绿色的点代表对照集合)，右方区域显示了迭代过程中的实时配准结果，右侧点云的外部轮廓表示点云的最小有向包围盒。

上述迭代过程结束后，便得到了最终的配准结果。为了让实验结果更加直观，图4展示了配准结果的部分放大图。如图4所示，图像左侧区域显示了原始集合与对照集合的初始位置，能够观察到它们之间存在较大的错位。左侧图像显示了配准前连个数据集合重叠区域的位置关系，右侧图像显示了本文算法的配准结果。经过本发明算法处理过后的结果显示在图像的右侧，图像显示大部分存在错位的点都得到了较为精确的配准。

在配准过程结束后，我们将得到的一系列变换矩阵作用于原始集合与对照集合上，并将它们用imageJ渲染出来。图5可视化了原始数据以及经过配准后的结果数据。图5渲染了四个眼底光学相干性断层扫描图像体数据。上面两幅图像分别是上述提到的原始集合与对照集合，后两幅本发明的实验结果，左下图像为配准结果的侧视图，右下为配准结果的俯视图，图中凹陷部分为人眼黄斑中心。采用传统的眼底体数据配准算法往往导致配准后的视网膜内核层、感光细胞层和色素上皮层出现断崖的现象。然而从图5显示结果来看，配准结果并未出现明显的视网膜层断层的情况，且在数据的重叠区域也观察不到清晰的镶嵌痕迹，这种具有较大视野的光学相干性断层扫描体数据可能对临床医师预防和诊断眼科疾病有一定的帮助。

在本发明中，图3-5仅用来表示本发明所获得的效果，但本发明的效果不仅依赖于图3-5来表示。

实施例2

此外，本发明采用了配准误差来评估算法的配准精确度，使用时间消耗长度来评估算法的配准效率。配准误差表示配准过程中对应点匹配失败的点数占总点云点数的百分比，它具有以下计算形式

其中Success(P_A,Q_B)的定义为

上述公式中N代表所有参与计算的对应点对个数。(P_A,Q_B)代表一个对应点对。Success(P_A,Q_B)用来表示对应点对(P_A,Q_B)的配准结果。若配准后对应点对间的欧氏距离小于给定的阈值δ，则表示此对应点对配准成功，Success(P_A,Q_B)＝1。本发明算法中δ的取值为(P_A,Q_B)精确配准前欧氏距离的0.15倍，即配准后欧氏距离小于原来距离15％的对应点对被视为配准成功。通过统计与计算，表1显示了本发明的算法与Besl提出的传统迭代最近点算法在时间消耗和配准误差方面的对比。

表1本发明算法与传统迭代最近点算法的实验结果的对比

作为实验对比，本发明对一些开源的点云数据集进行了实验。图6展示了采用本发明算法对公开库“斯坦福兔子”点云数据的配准结果，左图为点云数据的初始位置，右图为配准后的点云数据。此点云数据共有9731个点云点，使用传统迭代最近点算法用时11.023秒，配准误差为0.001130。使用改进的迭代最近点算法用时3.994秒，配准误差为0.000794。本发明通过对传统迭代最近点算法的改进，使得点云数据的配准从计算时间和配准精度上都得到了较为明显的提高。

为了更直观的对比改进算法和传统算法，本发明对不同大小的点云数据进行了对比实验，并采用折线图的方式展示了算法在计算时间上的差异。表2中第一个点云数据量为9731的数据集合为“斯坦福兔子”，第二个为斯坦福三维扫描仓库中的“龙”(http://graphics.stanford.edu/data/3Dscanrep/)，其余为实际项目中的点云数据。下面的表2和图7显示了对于不同数据集本发明算法与传统迭代最近点算法在配准时间上的对比。

表2不同大小的点云数据集合，本发明算法与传统迭代最近点算法计算时间对比

如表2与图7所示，当处理具有较大数据量的密集型点云时，本发明算法在时间复杂度和配准精度方面有明显的优势。在大多数情况下，本发明算法将传统的迭代最近点算法效率提高了70％。

上述实验结果验证了本发明对OCT眼底图像配准与拼接的有效性，同时验证了本发明对创建具有较大视野眼底视网膜图像的精确性和鲁棒性。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施方式仅限于此，对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单的推演或替换，都应当视为属于本发明由所提交的权利要求书确定保护范围。

Claims (9)

1.一种OCT眼底图像数据配准方法，包括如下步骤：

步骤一、采用Canny边缘检测方法对图像进行去噪、边缘检测：

对原始图像进行灰度化，对图像进行滤波，然后计算图像的梯度幅值，对梯度幅值进行非极大值抑制，然后用双阈值算法检测和连接边缘；

步骤二、图像点云数据的提取及可视化：

将经过步骤一处理的图像按照原有的顺序叠加成一个三维体数据，并把每一个边缘点按照其所处的空间位置抽取成三维点云中的一个点，点云点的三维坐标对应叠加得到的三维体数据的空间坐标；

步骤三、点云数据边界特征点的提取：

将点云中的点按照其空间坐标分配到不同的空间网格中，接着找出所有属于点云边界的空间网格，最后将边界网格中的点云点提取出来作为点云数据的边界特征点；

步骤四、采用奇异值分解算法完成点云数据初始配准：

令P代表原始集合，Q代表对照集合，参见公式3，定义目标矩阵如下：

C_P和C_Q分别是原始集合P以及对照集合Q的质心，M代表点云数据集中点云的数量，P_i和Q_i分别代表原始集合P与对照集Q合中第i个点，将目标矩阵E采用奇异值分解(SingularValue Decomposition,SVD)，则E＝UDV^T。其中U的列是EE^T矩阵的特征向量，V的列是E^TE矩阵的特征向量。V^T是V的转置矩阵且D为对角矩阵，D＝diag(d_i)，其中，d_i为E的奇异值，令

则旋转矩阵R＝UBV^T，平移矩阵T＝C_Q-RC_P，将所求得的矩阵R和T作用于原始集合P以消除初始条件下原始集合P与对照集合Q可能存在的较大位移误差；

步骤五、采用改进迭代最近点算法完成点云数据精确配准：

步骤5.1初始化。对于给定的两个点云集合P和Q，指定一个收敛阈值τ，

步骤5.2采用公式5算出点云数据中每个点对应的权值，比较阈值ε并将权值较低的点排除，

Q_B为集合Q中点P_A的对应点，Dis(P_A,Q_B)代表P_A和Q_B间的欧氏距离，Dis_MAX代表点对间欧氏距离的最大值，欧式距离采用公式6计算得到，

步骤5.3迭代以下步骤，直到公式6的最小均方根误差收敛于所给定的阈值τ：

步骤5.3.1根据公式6计算出集合P和Q中点云之间的欧式距离，

其中，ω_x，ω_y，ω_z分别代表M-估计在各个坐标方向的权重因素，(x_A,y_A,z_A)，(x_B,y_B,z_B)分别为集合P中点A和Q中点B的空间坐标，

步骤5.3.2对于排除了权值较低的点的集合P，在集合Q中寻找欧氏距离最近的点作为对应点并存储在最近点集中，

步骤5.3.3利用公式7采用最小二乘法计算计算集合P和最近点集间的旋转矩阵R和平移矩阵T，

步骤5.3.4将旋转矩阵R和平移矩阵T应用于集合P，得到新的集合，利用公式7计算最小均方根误差是否收敛于所给定的阈值τ，如果是则结束运算，否则利用步骤5.3进行迭代计算。

2.根据权利要求1所述的OCT眼底图像数据配准方法，其特征在于：

在步骤一中，原始图像的灰度化公式为Gray＝0.299R+0.587G+0.114B，图像滤波采用高斯滤波，采用一阶有限差分来计算图像像素矩阵关于横向和纵向的偏导数。

3.根据权利要求2所述的OCT眼底图像数据配准方法，其特征在于：

在步骤一中，采用双阈值来检测并连接边缘，高阈值的作用在于抑制边缘，所有高于此阈值的梯度幅值才能被视为边缘，低阈值的作用在于连接边缘，将所有梯度幅值高于低阈值的点视为边缘并连接起来成为最终边缘检测的结果。

4.根据权利要求1所述的OCT眼底图像数据配准方法，其特征在于：

在步骤三中，采用有向包围盒作为点云数据的最小包围盒，以将点云数据点划分到不同的空间网格中，

按照相同的体积将点云数据分解成若干个空间网格，每个空间网格的大小被定义为：

公式中L代表点云数据中点云点的数量，V代表最小包围盒的体积。V/L为点云密度的倒数，代表点云数据中每一个点云点所占空间的平均大小。令空间网格的初始大小S_grid为的点云密度倒数的K倍，并将点云数据最小包围盒按照S_grid大小分成若干个空间网格。

5.根据权利要求4所述的OCT眼底图像数据配准方法，其特征在于：

使用边界种子网格算法找到所有边界空间网格，并提取这个边界空间网格中包含的点云点作为点云数据的边界特征点，将空间网格分成两类，空格和实格，不包含任何点云点的空间网格为空格，其他的空间网格为实格。以空间坐标(x,y,z)表示某一网格空间位置，用函数f表示此空间网格的类型，如果某个空间网格是实格，f(x,y,z)＝1，否则f(x,y,z)＝0，使用公式2来判定某个空间网格是否为的边界空间网格：

当U(x,y,z)≤1时，代表此空间网格的六邻域中上下、左右、前后的网格中至多有一对都是实格，那么此空间网格是一个边界网格。

6.根据权利要求4或5所述的OCT眼底图像数据配准方法，其特征在于：

在步骤三中，8≤K≤24。

7.根据权利要求1所述的OCT眼底图像数据配准方法，其特征在于：

在步骤五中，收敛阈值τ＝0.2，比较阈值ε为，0.2≤ε≤0.4。

8.根据权利要求7所述的OCT眼底图像数据配准方法，其特征在于：

在步骤五中，M-估计在各个坐标方向的权重公式为

$w=\left\{\begin{array}{cc}1& \left|v\right|\≤c\\ \frac{c}{\left|v\right|}& \left|v\right|>c\end{array}\right.$

v为各个坐标方向上标准化的残差值，c为常数。

9.根据权利要求8所述的OCT眼底图像数据配准方法，其特征在于：

c＝1.345。