CN101593351A - 基于距离变换和刚性变换参数估计的眼底图像配准方法 - Google Patents

Abstract

本发明是一种基于距离变换和刚性变换参数估计的眼底图像配准方法，所述方法主要步骤包括：(1)加载眼底图像。(2)提取视盘中心，估计图像的平移参数。(3)对视盘邻域内的像素点，计算梯度向量，进行血管分割以及血管分布概率特征的计算，通过最小化两个概率分布的相对熵(Kullback－Leibler Divergence)得出图像旋转参数的估计。(4)对第三步中分割出来的血管，计算它们的欧式距离变换。(5)进行图像的精确配准。本发明方法是一种快速、精确、鲁棒的自动眼底图像配准算法，在眼底图像配准方面有着重大的应用价值。

Description

基于距离变换和刚性变换参数估计的眼底图像配准方法

技术领域

本发明涉及图像处理，模式识别技术，特别涉及一种基于距离变换以及刚性变换参数估计的眼底图像自动配准技术。

背景技术

目前主流的自动眼底图像配准方法，主要有基于血管特征的配准，基于互信息的配准方法等。这些方法都还存在一些不足之处，基于血管特征的配准，对一些质量比较低的图像，血管特征难以准确提取，基于互信息的配准，配准需要的时间比较长，而且可能会陷入局部极值，无法得到精确配准结果。因此，很多医生还是采用手动配准，手动配准的成功率和精度都比较高，但是它最大的缺点就是加大了医生的负担，配准的速度非常慢。

发明内容

针对现有技术的缺陷，本发明的目的是提供一种快速，精确，鲁棒，具有重大临床使用价值的自动眼底图像配准方法。该方法是一种基于距离变换和刚性变换参数估计的自动眼底图像配准方法，包括以下五个顺序处理模块：眼底图像加载模块，用于加载眼底图像；视盘中心提取模块，用来提取眼底图像的视盘中心，估计两幅图像的平移参数；视盘血管分割以及分布特征计算模块，用来分割视盘周围的血管，并计算血管的方向分布特征；视盘血管欧式距离变换计算模块，用来计算视盘周围血管的欧式距离变换；眼底图像精确配准模块，根据之前计算出的血管欧式距离变换，实现眼底图像的精确配准。

所述眼底图像加载模块，应用计算机编程语言C++编写程序，读取眼底图像，将图像转化为二维数组，存储在计算机中，以方便后续模块进行处理。

所述视盘中心提取模块，具体实现是应用计算机编程语言C++，根据视盘特征，编写程序，提取视盘中心，根据视盘中心在图像中的坐标位置，估计两幅图像的平移参数。

所述视盘血管分割以及分布特征计算模块，具体实现是应用计算机编程语言C++编写程序，实现如下功能：对视盘邻域内的像素点，利用2个5×5的方向掩模，计算梯度向量，进行血管分割以及血管分布概率特征的计算，通过最小化2个概率分布的相对熵(Kullback-Leibler Divergence)得出图像旋转参数的估计。

所述视盘血管欧式距离变换计算模块，具体实现是应用计算机编程语言C++编写程序，实现如下功能：对已经分割出来的视盘血管，计算他们的欧式距离变换，以实现图像的快速精确配准。

所述眼底图像精确配准模块，具体实现是应用计算机编程语言C++编写程序，实现如下功能：计算基于视盘血管欧式距离变换的图像相似性度量函数，利用梯度下降优化算法，寻找这个相似性度量函数的最大值，实现眼底图像的精确配准。

具体地，本发明提出一种基于距离变换和刚性变换参数估计的眼底图像配准方法，包括以下步骤：

1)加载第一眼底图像和第二眼底图像；

2)提取所述两幅眼底图像的视盘中心，估计所述两幅眼底图像的平移参数；

3)对所述两幅眼底图像进行视盘血管分割以及分布特征计算，估计所述两幅图像的旋转参数；

4)以所述第一眼底图像作为参考图像，对步骤3)中分割出来的血管像素点作为标志点集，计算以视盘中心为圆心，图像宽度的十分之一为半径的圆形邻域内的距离变换；

5)以所述第二眼底图像作为浮动图像，将步骤2)和3)中求得的刚性变换参数作为初始变换参数，将浮动图像中分割出来的血管像素点映射到基准图像上，计算两幅图像基于距离变换的相似性度量函数值，通过优化这个函数值，求得最佳的变换参数，完成两幅图像的精确配准。

进一步，步骤2)中所述视盘中心提取包括：

1)将图像的灰度最大值的85％设定为阈值；

2)对图像中灰度值高于阈值的像素点进行聚类，设定类内距离为图像宽度的十二分之一，同时设定类内像素点数阈值，得出候选的视盘中心点集；

3)对视盘中心点集中的候选点，分别计算其半径为图像宽度十二分之一的圆形区域内的血管密度，选择血管密度最大的点作为视盘中心。

进一步，步骤3)中所述进行视盘血管分割以及分布特征计算，估计两幅图像的旋转参数包括：

1)对以视盘中心为圆心，图像宽度的十二分之一为半径的圆形邻域内的像素点，利用两个5×5的方向掩模，计算梯度向量；

2)根据梯度向量进行血管分割以及血管方向概率分布特征的计算；

3)通过最小化两个血管方向概率分布的相对熵(Kullback-LeiblerDivergence)得出图像的旋转参数估计。

进一步，所述血管方向概率分布特征的计算步骤包括：对于分割出的血管边缘，计算梯度向量(G_x，G_y)，其中血管方向θ由如下表达式决定：

将整个圆周角[-π/2，3π/2]等分为72个方向，求出血管方向概率分布。

进一步，所述图像的旋转参数估计包括：对于已经求出的血管方向概率分布P₁，P₂，固定P₁，对P₂对应血管方向分布每次顺时针旋转一个方向，旋转71次，得出新的概率分布P₂ ¹，P₂ ²，...，P₂ ⁷¹，计算P₁与P₂，P₂ ¹，P₂ ²，...，P₂ ⁷¹的相对熵(Kullback-Leibler Divergence)，求出相对熵最小的上标i：

$i=\arg \underset{i}{\min }\underset{j=0}{\overset{71}{\mathrm{\Σ}}}{P}_{1_j}\log \frac{P_{1_j}}{{P_{2_j}}^i},$ i＝0，1，...71

${P_{2_j}}^0=P_{2_j},$ j＝0，1，...，71

下标j表示概率分布P在j方向上的分量，图像的旋转参数Δθ＝iπ/36。

进一步，步骤5)中所述两幅图像基于距离变换的相似性度量函数为：

$f=\frac{\underset{i=1}{\overset{\mathrm{numi}}{\mathrm{\Σ}}}\frac{1}{1+d_i^2}}{\mathrm{numi}}$

其中，numi是浮动图像分割出的血管边缘像素数目，d_i是浮动图像第i个血管边缘像素通过变换映射到参考图像空间上对应像素的欧式距离值，通过最大化函数值f，求取精确的图像变换参数，实现精确配准。

本发明利用图像特征预先估计图像的刚性变换参数，完成图像的粗配准，然后通过计算欧式距离变换，完成的图像的精确配准。由于预先准确估计图像变换参数，减少了大量的计算步骤，提高了算法的运行速度。同时，基于欧式距离变换的相似性度量函数对低质量的眼底图像也能很好的描述图像的匹配程度，算法具有比较高的鲁棒性。实验结果表明，本算法可以在1秒以内完成图像的配准，刚性配准精度可以达到亚像素级，对于多种质量低下的眼底图像可以完成自动配准。因此，具有重大的应用价值。

附图说明

图1是本实现本发明的模块示意图；

图2(a)是病人在2005年6月20日的眼底图像；

图2(b)是病人在2006年3月7日的眼底图像；

图2(c)是提取的图2(a)图视盘中心区域；

图2(d)是提取的图2(b)图视盘中心区域：

图2(e)是分割出的视盘邻域内的图2(a)图血管；

图2(F)是分割出的视盘邻域内的图2(b)图血管；

图3是血管方向划分示意图；

图4是欧式距离变换示意图，其中，像素值为0的点是标志点，欧式距离变换计算的是像素到最近的标志点的欧式距离；

图5(a)是2005年6月20日的图片作为参考图像；

图5(b)是2006年3月7日的图像作为浮动图像；

图5(c)是浮动图像经配准之后的结果图像；

图6是将参考图像与配准图像融合之后的新图像。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明的配准方法主要包括以下步骤：

步骤1：运行眼底图像加载模块101，加载眼底图像，应用计算机编程语言C++编写程序，读取眼底图像，将图像转化为二维数组，存储在计算机中，以方便后续模块进行处理。

步骤2：运行视盘中心提取模块102，估计两幅图像的平移参数。

视盘中心的提取依赖于视盘的以下特征属性：1)视盘一般对应着图像中最明亮的区域。2)视盘内血管众多，明显，血管密度是图像中最大的区域。首先，设定眼底图像中灰度最大值的85％为阈值T，设定类内距离L为图像宽度w的十二分之一，同时设定类内像素点数阈值N，对图像中灰度值高于阈值T的像素点进行聚类。初始类的数目设为0，类中心均设置为空。具体的聚类过程如下：

1、按顺序遍历图像的每一个高于阈值T的像素点，若当前类的数目为0，则将其作为第1个聚类中心，否则，计算当前像素点p与n个聚类中心的欧式距离L₁，L₂，...，L_n，选取其中距离最小值L_i，若L_i＜L，则将p归于第i类，更新第i类的中心坐标。否则将p作为第n+1类的聚类中心。

2、遍历完成后，统计每一类的元素数目n，只保留类内元素数目n大于N的类中心坐标，插入候选的视盘中心点集P。聚类完成，得到候选视盘中心点集P。

然后，对于点集P内的每一个候选点p_i，计算以p_i为中心，以图像宽度w的十二分之一为半径的圆形区域内的血管密度ρ。

血管密度ρ的计算步骤如下：

1、对于区域内的每一个像素点，利用2个5×5的gabor小波方向掩模：GradientX，GradientY，求解像素的梯度向量(G_x，G_y)。

$\mathrm{GradientX}=\left(\begin{array}{ccccc}0& -0.1338& 0& 0.1338& 0\\ 0& -0.1947& 0& 0.1947& 0\\ 0& -0.2206& 0& 0.2206& 0\\ 0& -0.1947& 0& 0.1947& 0\\ 0& -0.1338& 0& 0.1338& 0\end{array}\right)$

$\mathrm{GradientY}=\left(\begin{array}{ccccc}0& 0& 0& 0& 0\\ -0.1338& -0.1947& -0.2206& -0.1947& -0.1338\\ 0& 0& 0& 0& 0\\ 0.1338& 0.1947& 0.2206& 0.1947& 0.1338\\ 0& 0& 0& 0& 0\end{array}\right)$

2、对于梯度向量模大于经验阈值T_G的像素，认为它是血管的边缘，将其分割出来。

3、计算血管密度ρ， $\mathrm{\ρ}=\frac{\mathrm{num}}{S},$ 其中num是血管的边缘像素数，S是圆形区域的面积。

从这些候选点中，选取血管密度最大的点作为视盘中心，视盘中心区域为图2(c)，(d)白色部分，视盘中心在图上以黑色十字箭头标识，这样我们可以获得2幅眼底图像的视盘中心坐标，记为(x₁，y₁)，(x₂，y₂)，得出图像的平移参数估计：Δx＝x₁-x₂，Δy＝y₁-y₂。

步骤3：运行视盘血管分割及分布计算模块103，估计图像的旋转参数。

图像旋转参数的估计主要由以下几步完成：

1、按照步骤2的血管分割方法，分割出以视盘中心为圆心，图像宽度w的十二分之一为半径的圆形区域内的血管边缘像素，如图2(e)，(f)所示，然后计算每一点的血管方向θ。

$\mathrm{\&theta;}=\left\{\begin{array}{cc}{\tan }^{-1}\frac{G_y}{G_x}& \mathrm{if}{G}_x<0\\ {\tan }^{-1}\frac{G_y}{G_x}+\mathrm{\&pi;}& \mathrm{else}\end{array}\right.$

其中，(G_x，G_y)是由GradientX，GradientY这2个掩模运算求得的像素梯度向量。

2、将整个圆周角[-π/2，3π/2]等分为72个方向，如图3所示，求出血管方向概率分布P＝[P₀，P₁，...，P₇₁]。其中：numθ_i是第i个方向内的血管边缘像素数目。

$P_i=\frac{\mathrm{num}{\mathrm{\&theta;}}_i}{\underset{j=0}{\overset{71}{\mathrm{\&Sigma;}}}{\mathrm{num\&theta;}}_j},$ i＝0，1，...，71

3、在求出了两幅图像的血管方向概率分布P₁，P₂后，固定P₁，将P₂对应血管方向分布每次顺时针旋转一个方向，旋转71次，得出新的概率分布P₂ ¹，P₂ ²，…，P₂ ⁷¹，计算P₁与P₂，P₂ ¹，P₂ ²，…，P₂ ⁷¹的相对熵(Kullback-Leibler Divergence)，求出相对熵最小的上标i：

$i=\arg \underset{i}{\min }\underset{j=0}{\overset{71}{\mathrm{\&Sigma;}}}{P}_{1_j}\log \frac{P_{1_j}}{{P_{2_j}}^i},$ i＝0，1，...71

${P_{2_j}}^0=P_{2_j},$ j＝0，1，...，71

下标j表示概率分布P在j方向上的分量。则图像的旋转参数

Δθ＝iπ/36。

步骤4：以任意一幅图像作为参考图像，运行视盘血管欧式距离变换计算模块104。将步骤3中分割出来的血管像素点作为标志点集，计算以视盘中心为圆心，图像宽度的十分之一为半径的圆形邻域内的欧式距离变换。距离变换是图像处理中的常用技术，给定一幅图像及标志像素点集，如何求出图像中的每个像素点到最近的标志点的距离就是距离变换所要解决的问题。一幅已经完成欧式距离变换计算的示意图像如图4所示。图像大小是8×7，像素值为0的点对应的就是标志点，其余的像素值代表该像素距离最近标志点的欧式距离。

步骤5：将另一幅图像作为浮动图像，运行眼底图像精确配准模块105。将步骤2，3求得的刚性变换参数作为初始变换参数，将浮动图像中分割出来的血管像素点映射到基准图像上，计算两幅图像基于距离变换的相似性度量函数值f：

$f=\frac{\underset{i=1}{\overset{\mathrm{numi}}{\mathrm{\&Sigma;}}}\frac{1}{1+d_i^2}}{\mathrm{numi}}$

其中，numi是浮动图像分割出的血管边缘像素数目。d_i是浮动图像第i个血管边缘像素通过变换映射到参考图像空间上对应像素的欧式距离值。我们通过基于梯度下降的寻优算法，快速寻找相似性度量函数f的最大值，解得相应的变换参数，实现图像的快速，精确配准。梯度下降的寻优算法是一种常用的优化算法，主要的算法流程如下：

对需要优化的函数f(x)在x_k进行泰勒展开，省略二阶及以上小量：

$f(x_k+d)=f\left(x_k\right)+d^T\&dtri;f\left(x_k\right)+...$

其中d代表从x_k点开始的第k+1步迭代的寻优方向，在梯度下降算法中，d选取为x_k点的负梯度方向，第k+1步寻优的表达式如下：

$x_{k+1}=x_k+{\mathrm{\&alpha;}}_k(-\&dtri;f\left(x_k\right))$

α_k是第k+1步迭代的步长，当函数的梯度值趋于0时，停止迭代，优化过程结束，此时的x_k即为所求的最优解。

为了验证本发明方法，我们选取了临床获得的眼底图像作为实验样本，通过大量的实验，我们可以在1秒以内完成亚像素级的刚性眼底图像配准，对低质量的眼底图像也能得到很好的实验结果。具体的配准结果如图5(c)，图6所示。其中图6是将配准后的图像与之前的参考图像融合之后得到的新图像，图像中亮度高的矩形区域来源于配准后的图像，亮度低的矩形区域来源于参考图像。从图上可以看出，配准的精度非常高。实验表明，本发明是快速，精确，鲁棒的，具有巨大的实用价值。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (6)

1.一种基于距离变换和刚性变换参数估计的眼底图像配准方法，包括以下步骤：

1)加载第一眼底图像和第二眼底图像；

2)提取所述两幅眼底图像的视盘中心，估计所述两幅眼底图像的平移参数；

3)对所述两幅眼底图像进行视盘血管分割以及分布特征计算，估计所述两幅图像的旋转参数；

4)以所述第一眼底图像作为参考图像，对步骤3)中分割出来的血管像素点作为标志点集，计算以视盘中心为圆心，图像宽度的十分之一为半径的圆形邻域内的距离变换；

5)以所述第二眼底图像作为浮动图像，将步骤2)和3)中求得的刚性变换参数作为初始变换参数，将浮动图像中分割出来的血管像素点映射到基准图像上，计算两幅图像基于距离变换的相似性度量函数值，通过优化这个函数值，求得最佳的变换参数，完成两幅图像的精确配准。

2.根据权利要求1所述的方法，其特征在于，步骤2)中所述视盘中心提取包括：

1)将图像的灰度最大值的85％设定为阈值；

2)对图像中灰度值高于阈值的像素点进行聚类，设定类内距离为图像宽度的十二分之一，同时设定类内像素点数阈值，得出候选的视盘中心点集；

3)对视盘中心点集中的候选点，分别计算其半径为图像宽度十二分之一的圆形区域内的血管密度，选择血管密度最大的点作为视盘中心。

3.根据权利要求1所述的方法，其特征在于，步骤3)中所述进行视盘血管分割以及分布特征计算，估计两幅图像的旋转参数包括：

1)对以视盘中心为圆心，图像宽度的十二分之一为半径的圆形邻域内的像素点，利用两个5×5的方向掩模，计算梯度向量；

2)根据梯度向量进行血管分割以及血管方向概率分布特征的计算；

3)通过最小化两个血管方向概率分布的相对熵(Kullback-LeiblerDivergence)得出图像的旋转参数估计。

4.根据权利要求3所述的方法，其特征在于，所述血管方向概率分布特征的计算步骤包括：对于分割出的血管边缘，计算梯度向量(G_x，G_y)，其中血管方向θ由如下表达式决定：

将整个圆周角[-π/2，3π/2]等分为72个方向，求出血管方向概率分布。

5.根据权利要求3所述的方法，其特征在于，所述图像的旋转参数估计包括：对于已经求出的血管方向概率分布P₁，P₂，固定P₁，对P₂对应血管方向分布每次顺时针旋转一个方向，旋转71次，得出新的概率分布P₂ ¹，P₂ ²，...，P₂ ⁷¹，计算P₁与P₂，P₂ ¹，P₂ ²，...，P₂ ⁷¹的相对熵(Kullback-LeiblerDivergence)，求出相对熵最小的上标i：

$i=\arg \underset{i}{\min }\underset{j=0}{\overset{71}{\mathrm{\&Sigma;}}}{P}_{1_j}\log \frac{P_{1_j}}{{P_{2_j}}^i},$ i＝0，1，...71

${P_{2_j}}^0=P_{2_j},$ j＝0，1，...，71

下标j表示概率分布P在j方向上的分量，图像的旋转参数Δθ＝iπ/36。

6.根据权利要求1所述的方法，其特征在于，步骤5)中所述两幅图像基于距离变换的相似性度量函数为：

$f=\frac{\underset{i=1}{\overset{\mathrm{numi}}{\mathrm{\&Sigma;}}}\frac{1}{1+d_i^2}}{\mathrm{numi}}$

其中，numi是浮动图像分割出的血管边缘像素数目，d_i是浮动图像第i个血管边缘像素通过变换映射到参考图像空间上对应像素的欧式距离值，通过最大化函数值f，求取精确的图像变换参数，实现精确配准。

Images