基于图像配准的多模态视网膜图像融合方法及系统与流程

本发明涉及图像处理领域，具体而言，涉及基于图像配准的多模态视网膜图像融合方法及系统。

背景技术：

眼底视网膜图像是诊断包括青光眼、年龄相关性黄斑变性在内的多种视网膜疾病的重要依据。此外，眼底作为人体唯一一个可直视的血管窗口，可在一定程度上反映全身其他脏器的血流动力学改变情况。多模态视网膜图像通常包含重要的视网膜的局部结构和时相信息。如果能将同一患者的不同模态的视网膜图像进行融合，可以为医生提供病变部位的多种互补信息，从而为医生决策诊疗提供更为全面清晰的依据。

由于眼底图像血管网络结构复杂，其融合效果受到低质量视网膜图像的影响，包括不均匀的内容对比度和非线性强度差异，以及非血管、无纹理区域和各种病理病变区域导致的重叠恶化。因此，高精度的视网膜图像配准是解决上述问题的关键步骤。现有的视网膜图像配准方法大致分为两类：基于区域的配准方法和基于特征的配准方法。基于特征的方法通常包括提取特征和构建变换两个步骤，以估计图像特征之间的对应关系。

最经典的基于特征的配准方法是采用软分配和确定性退火分别估计对应关系和控制薄板样条(tps)转换更新的tps-rpm方法。此后，基于经典的tps-rpm方法，有文献提出了一种基于全局和局部混合距离(glmdtps)的鲁棒方法用于点集配准问题。

虽然上述方法在特征匹配、点集和图像配准的应用方面取得了一些成功，但仍然存在一些问题。首先，由于视网膜图像的特殊性(重复结构、多模态)，不准确的特征提取和描述会产生大量错误的特征匹配。此外，在图像空间变换过程中，上述研究没有采用约束或仅采用单一约束的方法会导致图像空间变换的不准确。这些问题都将影响基于图像配准的多模态视网膜图像融合的稳健性和精确性。

综上，现有的视网膜图像配准方法普遍存在特征提取和描述不准确导致产生大量错误的特征匹配，以及约束方法导致图像控件变换的不准确，进而影响多模态视网膜图像融合的稳健性和精确性。因此，需要一种多模态视网膜图像融合方法，能够解决上述问题。

技术实现要素：

基于现有技术存在的问题，本发明提供了一种基于图像配准的多模态视网膜图像融合方法及系统。具体方案如下：

一种基于图像配准的多模态视网膜图像融合方法，包括：

图像输入：获取包括源图像和目标图像的视网膜图像对；

图像预处理：对所述视网膜图像对进行预处理，获取视网膜边缘图像对，从每对所述视网膜边缘图像对中提取特征点集，所述特征点集包括从预处理后的源图像中提取的源特征点集和从预处理后的目标图像中提取的目标特征点集；

特征组合和提取：组合多个特征描述符构建多特征差异描述子，引导对所述特征点集的特征提取；

特征点集配准：通过所述多特征差异描述子评估所述源特征点集与所述目标特征点集之间的对应关系，获取对应关系矩阵，根据所述对应关系矩阵更新所述源特征点集的空间变换；

图像配准：根据空间变换前后的源特征点集实现图像配准，结合所述目标图像获取空间变换后的图像；

图像融合：将空间变换后的图像和预设参考图像进行像素级的融合得到融合图像。

在一个具体实施例中，多个所述特征描述符包括基于类sift尺度空间的边缘定向直方图描述符、局部几何结构特征描述符、全局几何结构特征描述符；

所述特征组合和提取具体包括：

通过所述基于类sift尺度空间的边缘定向直方图，匹配位于相同场景、不同光谱条件下的所述视网膜边缘图像对上的特征点集，获取所述源特征点集和所述目标特征点集的尺度差异；通过所述局部几何结构特征描述符，获取所述源特征点集和所述目标特征点集的局部特征差异矩阵；通过所述全局几何结构特征描述符，获取所述源特征点集和所述目标特征点集的全局特征差异矩阵；根据所述尺度差异、所述局部特征差异矩阵和所述全局特征差异矩阵构建所述多特征差异描述子。

在一个具体实施例中，“通过所述基于类sift尺度空间的边缘定向直方图，匹配位于相同场景、不同光谱条件下的所述视网膜边缘图像对上的特征点集，获取所述源特征点集和所述目标特征点集的尺度差异”具体包括：

通过基于类sift尺度空间表示来检测所述特征点集；通过边缘定向直方图描述符描述所述特征点集，获取所述特征点集的特征描述，所述边缘定向直方图包含来自每个特征点附近的轮廓的空间信息，以每个特征点为中心设置预设尺寸的窗口描述所述视网膜边缘图像对的形状和轮廓，并维护尺度空间的不变性；根据所述特征描述，匹配所述目标特征点集中的特征点与所述源特征点集中的特征点，通过使用欧几里得距离来评判特征点之间的差异程度，获取所述尺度差异。在一个具体实施例中，所述基于类sift尺度空间的边缘定向直方图描述符通过尺度差异来提高特征描述的鲁棒性；

所述尺度差异的定义为：

sd(x，y)＝sclx-scly,(x＜sd＜y)

其中，sd表示所述尺度差异，scl是特征点所处于的尺度空间，x表示源特征点集，y表示目标特征点集，x和y的定义如下：

其中，表示sd直方图中的峰值。

在一个具体实施例中，假设特征点集t包含j个特征点，则定义tik是特征点ti的第k个最近相邻特征点，则点集t的第i个点ti的局部几何特征描述符定义如下：

其中，k表示相邻点的数量，vik表示从ti到tik的向量，μik是一个权重参数，用于控制相邻点矢量对局部特征的描述力度；μik的定义如下：

所述局部特征差异矩阵中的每个元素定义为：

其中，ψlmn表示所述局部特征差异矩阵中某一元素，m表示目标特征点集中的特征点，n表示源目标点集中的特征点，表示高斯径向基函数，y表示目标特征点集。

在一个具体实施例中，通过欧几里得距离描述所述全局特征差异，所述全局特征差异矩阵中的每个元素定义为：

ψgmn＝||ym-c(ρn，θn)||

其中，ψgmn表示所述全局特征差异矩阵中某一元素，m表示目标特征点集中的特征点，n表示源特征点集中的特征点，c(ρn，θn)表示空间变换后高斯混合模型的质心。

在一个具体实施例中，所述多特征差异描述子定义如下：

mlf(t)＝ψg+sd(x,y)+t1ψl

其中，t1是所述局部特征描述符的退火参数，ψg表示所述全局特征差异矩阵，ψl表示所述局部特征差异矩阵，sd(x,y)表示尺度差异，mlf(t)表示所述多特征差异描述子。

在一个具体实施例中，所述特征点集配准具体包括：

在所述源特征点集上构造概率模型，通过所述多特征差异描述子计算所述源特征点集与所述目标特征点集之间的对应关系矩阵；基于所述对应关系矩阵，调整所述概率模型的参数，基于全局几何结构约束和局部几何结构约束更新所述源特征点集的空间变化；迭代上述步骤，使所述源特征点集逐步接近所述目标特征点集，直至满足预设匹配关系。

在一个具体实施例中，“在所述源特征点集上构造概率模型，通过所述多特征差异描述子计算所述源特征点集与所述目标特征点集之间的对应关系矩阵”具体包括：

在所述源特征点集上构造高斯混合模型；通过所述高斯混合模型测量所述源特征点集与所述目标特征点集之间的所述多特征差异描述子来评估对应关系；将所述对应关系的评估问题转换为所述高斯混合模型的概率密度，使用近似解求解所述概率密度；基于所述概率密度，根据贝叶斯规则计算所述高斯混合模型的后验概率来估计所述对应关系，获取所述对应关系矩阵；根据所述对应关系矩阵和更新前的源特征点集获取更新后的源特征点集。

在一个具体实施例中，“基于所述对应关系矩阵，调整所述概率模型的参数，基于全局几何结构约束和局部几何结构约束更新空间变化”具体包括：

基于所述对应关系矩阵，通过调整所述高斯混合模型的模型参数更新空间变化，获取所述模型参数的最优解；添加所述全局几何结构约束，维持所述源特征点集在空间变化时的全局稳定性；基于所述局部几何特征描述符，添加所述局部几何结构约束，通过判断空间变换前后的特征点集之间的局部相似性，约束所述源特征点集的局部变形。

在一个具体实施例中，所述模型参数包括第一优化参数σ²和第二优化参数θ；所述模型参数的最优解的获取过程包括：

通过最小化所述概率密度的负对数似然函数获取所述模型参数的最优解；采用期望最大化算法求解所述概率密度的负对数似然函数，所述期望最大化算法包括计算所述后验概率和最大化期望；根据所述后验概率，将所述概率密度的负对数似然函数作为所述最大化期望的能量函数；基于所述全局几何结构约束，在所述能量函数中加入全局约束项，获取第一约束函数；基于所述局部几何结构约束，在所述第一约束函数中加入局部约束项，进而获取所述高斯混合模型的参数的最优解。

在一个具体实施例中，所述概率密度的表达式为：

表示高斯混合模型的概率密度函数，pdf(ym|ρn)表示一个高斯组件ρn中对数据点ym的概率密度函数，其中包括n+1个高斯组件，m表示目标特征点集中的特征点，n表示源目标点集中的特征点，m表示目标特征点集中的特征点总数，n表示源目标点集中的特征点总数，每个高斯组件ρn表示源特征点集x中的xn，p(n)＝1/a，额外的第n+1个高斯组件用于消除冗余点的影响，使冗余点受参数ω约束(0＜ω＜1)；

一个高斯组件ρn中对数据点ym的概率密度如下：

pdf(ym|ρn)表示一个高斯组件ρn中对数据点ym的概率密度，mlf(t)表示尺度差异，σ²表示第一优化参数；所述对应关系矩阵表达式如下：

p^old(ρn|ym)表示一个高斯组件ρn中对数据点ym的对应关系矩阵，n表示源目标点集中的特征点，n表示源目标点集中的特征点总数。

在一个具体实施例中，所述概率密度的负对数似然函数表达式为：

e表示概率密度的负对数似然函数，ρn表示一个高斯组件中，ym为数据点，c(ρn，θn)表示空间变换后高斯混合模型的质心，p(x)表示对应关系；

将所述概率密度的负对数似然函数作为所述最大化期望的能量函数，展开如下：

p^old(ρn|ym)表示一个高斯组件ρn中对数据点ym的对应关系矩阵，θ表示第二优化参数；基于高斯径向基函数定义源特征点集x的空间变换：

其中，表示源特征点集的空间变换，w是高斯内核的n×d权重矩阵，g是由高斯径向基函数获得的n×n高斯内核矩阵，该矩阵中的某一元素的表达式为：

其中，gnm为高斯内核矩阵中的元素，n、m表示该元素在矩阵中的位置，x表示源特征点集中的特征点；其中，w是高斯内核的n×d权重矩阵，将第二优化参数θ转化为w，能量函数的表达式如下：

q(w,σ²)表示能量函数，p^old(ρn|ym)表示一个高斯组件ρn中对数据点ym的对应关系矩阵。

在一个具体实施例中，所述全局几何结构约束包括在所述能量函数中添加全局几何约束项，所述全局几何约束项为一个正则化运算符，所述正则化运算符的表达式为：

其中，常量λ是控制全局约束强度的权重参数，r表示正则化运算符，表示源特征点集的空间变换，trace()表示矩阵的迹，w是高斯内核的n×d权重矩阵，g是由高斯径向基函数获得的n×n高斯内核矩阵；

在所述能量函数中添加全局几何约束项获取第一能量函数，表达式为：

qg(w,σ²)表示第一能量函数，q(w,σ²)表示能量函数，常量λ是控制全局约束强度的权重参数。

在一个具体实施例中，所述局部几何结构约束项表达式如下：

其中，η是控制局部约束强度的权重参数，通过在空间变换之前和之后判断x的局部相似性来约束空间变换中源特征集x的局部变形；通过确定性退火技术定义η，表达式如下：

其中，m是预先设置的权重参数η的最大值，t是当前的迭代次数，而常量c用于控制约束强度变化的速度。

在一个具体实施例中，在所述第一能量函数表达式中添加局部几何结构约束项获取第二能量函数，表达式为：

qgl(w,σ²)＝qg(w,σ²)+ηtrace(d^td)

由此得到新的公式：

由于每一项的推导：

所述模型参数的最优解为：

w＝(gpxg+2ησ²g^tu^tug)^-1(gp0y-gp0y)

px为n×n矩阵，py为m×m矩阵，px和py是对角矩阵，分别由列向量p01和组成，1是一个元素都是1的列向量。

在一个具体实施例中，“根据空间变换前后的源特征点集实现图像配准，结合所述目标图像获取空间变换后的图像”具体包括：

根据所述对应关系矩阵获取空间变换后的源特征点集；根据空间变换后的源特征点集和空间变换前的源特征点集构建控制点集，实现图像配准；基于反推法，以所述目标图像和所述控制点集作为控制点，根据空间变换反推出变换后的图像。

一种基于图像配准的多模态视网膜图像融合系统，包括：

图像输入单元：用于获取包括源图像和目标图像的视网膜图像对；

图像预处理单元：用于对所述视网膜图像对进行预处理，获取视网膜边缘图像对，从每对所述视网膜边缘图像对中提取特征点集，所述特征点集包括从预处理后的源图像中提取的源特征点集和从预处理后的目标图像中提取的目标特征点集；

特征组合和提取单元：用于组合多个特征描述符构建多特征差异描述子，引导对所述特征点集的特征提取；

特征点集配准单元：用于通过所述多特征差异描述子评估所述源特征点集与所述目标特征点集之间的对应关系，获取对应关系矩阵，根据所述对应关系矩阵更新所述源特征点集的空间变换；

图像配准单元：用于根据空间变换前后的源特征点集实现图像配准，结合所述目标图像获取空间变换后的图像；

图像融合单元：用于将空间变换后的图像和预设参考图像进行像素级的融合得到融合图像。

在一个具体实施例中，所述特征点集配准单元具体包括，

模型构建单元：用于在所述源特征点集上构造高斯混合模型，通过所述多特征差异描述子计算所述源特征点集与所述目标特征点集之间的对应关系矩阵；

参数调整单元：用于基于所述对应关系矩阵，调整所述高斯混合模型的参数，基于全局几何结构约束和局部几何结构约束更新所述源特征点集的空间变化；

迭代单元：用于迭代所述模型构建单元和所述参数调整单元，使所述源特征点集逐步接近所述目标特征点集，直至满足预设匹配关系。

在一个具体实施例中，多个所述特征描述符包括基于类sift尺度空间的边缘定向直方图描述符(eoh-sift描述符)、局部几何结构特征描述符、全局几何结构特征描述符；

所述特征组合和提取单元具体包括：

尺度差异单元：用于通过所述基于类sift尺度空间的边缘定向直方图，匹配位于相同场景、不同光谱条件下的所述视网膜边缘图像对上的特征点集，获取所述源特征点集和所述目标特征点集的尺度差异；

局部特征单元：用于通过所述局部几何结构特征描述符，获取所述源特征点集和所述目标特征点集的局部特征差异矩阵；

全局特征单元：用于通过所述全局几何结构特征描述符，获取所述源特征点集和所述目标特征点集的全局特征差异矩阵；

差异构建单元：用于根据所述尺度差异、所述局部特征差异矩阵和所述全局特征差异矩阵构建所述多特征差异描述子。

本发明具有如下有益效果：

本发明针对现有技术中视网膜图像配准时匹配精度不足和空间变换不准确的缺陷，提出了一种基于图像配准的多模态视网膜图像融合方法及系统。基于视网膜图像的特殊性，采用基于特征的配准方法，以估计图像特征之间的对应关系，实现准确的特征提取和特征描述。在图像空间变换过程中，采用局部约束和全局约束相结合的约束方法，保证图像变换过程中的稳定性和准确性。本发明方法能够获得最佳的配准性能，并且在大多数情况下都优于目前最先进的方法。在高精度图像配准的基础上，进行有效的多模态眼底图像融合，便于眼科医生观察病变区域的变化，进一步协助相关视网膜疾病的诊断和治疗。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1是本发明实施例1的多模态视网膜图像融合方法流程图；

图2是本发明实施例1的输入数据示意图；

图3是本发明实施例1的视网膜图像融合流程示意图；

图4是本发明实施例1的视网膜图像融合流程示意图；

图5是本发明实施例1的视网膜图像融合流程示意图；

图6是本发明实施例1的多模态视网膜图像融合方法原理图；

图7是本发明实施例2的多模态视网膜图像融合系统示意图。

附图标记：

1-图像输入单元；2-图像预处理单元；3-特征组合和提取单元；4-特征点集配准单元；5-图像配准单元；6-图像融合单元；31-尺度差异单元；32-局部特征单元；33-全局特征单元；34-差异构建单元；41-模型构建单元；42-参数调整单元；43-迭代单元。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

基于图像配准的多模态视网膜图像融合方法具体包括：首先，使用引导图像滤波来增强不同类型的视网膜图像的边缘。其次，采用多特征描述符来组合视网膜图像的边缘图和几何结构特征，以改进特征集的描述和排除冗余点。然后，多特征引导模型为特征集配准提供准确的引导，并且基于特征点的图像配准最终提供准确的图像配准。最后，基于准确的图像配准进行多模态的视网膜眼底图像融合。

实施例1

本实施例提出了一种基于图像配准的多模态视网膜图像融合方法，如说明书附图1-6所示。流程步骤如说明书附图1，具体方案如下:

s1、图像输入：获取包括源图像和目标图像的视网膜图像对；

s2、图像预处理：对视网膜图像对进行预处理，获取视网膜边缘图像对，从每对边缘图像对中提取特征点集，特征点集包括从预处理后的源图像中提取的源特征点集和从预处理后的目标图像中提取的目标特征点集；

s3、特征组合和提取：组合多个特征描述符构建多特征差异描述子，引导对特征点集的特征提取；

s4、特征点集配准：通过多特征差异描述子评估源特征点集与目标特征点集之间的对应关系，获取对应关系矩阵，根据对应关系矩阵更新源特征点集的空间变换；

s5、图像配准：根据空间变换前后的源特征点集实现图像配准，结合目标图像获取空间变换后的图像；

s6、图像融合：将空间变换后的图像和预设参考图像进行像素级的融合得到融合图像。

具体地，s1、获取视网膜图像对，每对视网膜图像对包括源图像和目标图像。视网膜图像对如说明书附图2所示。

s2、图像预处理：对视网膜图像对进行边缘处理，获取视网膜边缘图像对，在视网膜边缘图像对中提取特征点集，特征点集包括源特征点集和目标特征点集。具体地，对视网膜图像对进行边缘处理，获取视网膜边缘图像对，引导图像滤波来增强不同类型的视网膜图像的边缘。视网膜图像对中包括源图像和目标图像，边缘视网膜图像对中也包括源边缘图像和目标边缘图像。在源边缘图像中提取特征点，获取源特征点集，在目标边缘图像中提取特征点，获取目标特征点集。将源特征点集记为将目标特征点集记为本实施例通过对两个特征点集进行配准来间接完成图像配准。

在视网膜图像配准中使用边缘图有三个主要原因和优势。1.通常整个图像中都分布有边缘。2.视网膜图像的边缘信息稳定且易于检索。3.边缘图中忽略了视网膜图像强度的梯度方向，由于其在多模态图像中的不变性，边缘图仅保留了图像的梯度幅度。通过采用边缘定向eoh-sift特征变换算法来可靠地提取独特且高度可重复的特征，以增强对比度并消除图像的噪声。

图像预处理的步骤如下：1.将强度直方图均衡为一个高斯分布，其中平均值mo＝124和方差so＝58，并使用引导滤波对所得图像进行去噪。2.sobel滤波器用于计算对比度并增强图像的边缘响应。3.使用对比度受限的自适应直方图均衡再次增强边缘的对比度。图像预处理方法处理视网膜图像可以产生大量的真实匹配，以用于对比度增强和噪声消除。

s3、特征组合和提取：将三个特征描述符组合到多特征差异描述子(mlf)引导特征提取。采用多特征描述符来组合视网膜图像的边缘图和几何结构特征，以改进特征集的描述和排除冗余点。三个特征描述符包括基于类sift尺度空间的边缘定向直方图(eoh-sift)描述符、全局几何结构特征和局部几何结构特征(lgsf)。三个特征描述符结合多特征差异描述子，实现互补，以增强描述符的性能并提高每个特征点的可识别性。

与sift相比，eoh-sift描述符使用尺度差异(sd)来提高特征描述的鲁棒性，通过丢弃在边缘图像中的子区域没有信息的特征点，即特征点所在的子区域中仅包含少量轮廓。给定一对特征集x和y，sd定义如下：

sd(x，y)＝sclx-scly,(x＜sd＜y)

其中，sd表示尺度差异，scl是特征点所处于的尺度空间，x表示源特征点集，y表示目标特征点集，x和y的粗略值通过如下两步计算获取：(1)计算所有匹配的sd的直方图；(2)提取sd直方图中的峰值

x和y的定义如下：

进一步的，所述的局部几何结构特征描述符(lgsf)是指：假设特征点集t包含j个特征点，则定义tik是特征点ti的第k个最近相邻特征点，则设置点集t的第i个点ti的局部特征描述符定义如下：

其中k表示相邻点的数量，vik表示从ti到tik的向量，μik是一个权重参数，用于控制相邻点矢量对局部特征的描述力度。由于向量特征包含欧几里德距离和方向的空间信息，本实施例通过采用和向量以全面的描述ti的局部特征。

局部几何结构特征描述符(lgsf)的性能主要受两个因素的影响：相邻点的数量k和权重参数μik。由于一个特征点ti的局部特征主要受其周围相邻点的影响，因此k的值应考虑到点周围可能产生影响的相邻特征点。例如在二维点集配准的情况下，应该至少考虑四个方向(上，下，左，右)以使lgsf具有较好的性能，即k的值应该取大于4的值。局部几何结构特征描述符(lgsf)可以应用于任何大于或等于二维的点集配准，因此k的值应该取决于具体应用的维度。在本实施例中，使用lgsf完成视网膜图像配准，即应用于二维情况，因此k取值为5。

向量vik的长度过长会对局部几何结构特征描述符(lgsf)的描述性能的产生较大负面影响，所以本实施例采用欧几里德范数定义每个向量的权重μik，并满足以下条件：矢量长度越短，影响越大。μik的定义如下：

lgsf的核心思想是特征点的局部特征可以由特征点与指定数量的相邻特征点之间的矢量之和来表示，并且这些矢量对lgsf描述子的贡献力度受其自身长度的影响。局部特征描述符对应的特征差异矩阵ψl的每个元素都可基于如下公式获取：

其中，ψlmn表示局部特征描述符对应的特征差异矩阵中的某个元素，x表示源特征点集，y表示目标特征点集，m表示目标特征点集中的特征点，n表示源目标点集中的特征点，表示高斯径向基函数。

同理，全局特征描述符对应的特征差异矩阵ψg的每个元素都可基于如下公式获取：

ψgmn＝||ym-c(ρn，θn)||

其中，ψgmn表示全局特征描述符对应的特征差异矩阵中的某个元素，c(ρn，θn)是空间变换后的高斯混合模型的质心，混合模型在s3进行了详细说明。

在本实施例中，将三个特征描述符组合到多特征差异描述子(mlf)中，mlf的表达式为：

mlf(t)＝ψg+sd(x,y)+t1ψl

其中，t1是局部特征描述符的退火参数，ψg表示全局特征描述符对应的特征差异矩阵，ψl表示局部特征描述符对应的特征差异矩阵，sd(x,y)表示尺度差异。

从一对视网膜图像中提取两个特征点集之后，通常存在大量冗余点，为了解决冗余点问题造成的配准问题，本实施例提出了一个鲁棒的多特征描述子引导的特征点集配准模型。给定两个特征点集，分别是从待配准图像提取的源特征点集以及从目标图像中提取目标特征点集特征点集配准模型包含以下两个主要步骤：(i)对应关系评估：在每次迭代时通过使用上述提出的mlf描述x和y以评估它们之间的对应关系；(ii)空间变换更新：根据(1)中评估的对应关系来构建x的非刚性转换以更新源点集的位置。通过迭代步骤(i)和(ii)使x逐渐接近目标y以匹配两点集之间对应的特征点。

s4、特征点集配准：在源特征点集上构造高斯混合模型，通过多特征差异描述子评估源特征点集与目标特征点集之间的对应关系，获取对应关系矩阵，根据对应关系矩阵更新源特征点集的空间变换。

具体地，特征点集配准具体包括：在源特征点集上构造高斯混合模型，通过多特征差异描述子计算源特征点集与目标特征点集之间的对应关系矩阵；基于对应关系矩阵，调整高斯混合模型的参数，基于全局几何结构约束和局部几何结构约束更新空间变化；迭代上述步骤，使特征点集逐步接近目标特征点集，直至满足预设匹配关系。

其中，在每次迭代时，通过多特征差异描述子描述源特征点集和目标特征点集之间的对应关系。每次空间变换后，源特征点集都会发生变化，逐渐接近目标特征点集。本实施例通过mlf描述变化的源特征点集和目标特征点集之间的对应关系。

高斯混合模型(gmm)就是用高斯概率密度函数(正态分布曲线)精确地量化事物，将一个事物分解为若干的基于高斯概率密度函数(正态分布曲线)形成的模型。对图像背景建立高斯模型的原理及过程：图像灰度直方图反映的是图像中某个灰度值出现的频次，也可以以为是图像灰度概率密度的估计。如果图像所包含的目标区域和背景区域相差比较大，且背景区域和目标区域在灰度上有一定的差异，那么该图像的灰度直方图呈现双峰-谷形状，其中一个峰对应于目标，另一个峰对应于背景的中心灰度。

对应关系评估是指首先使用高斯混合模型(gmm)通过测量两个特征点集之间的mlf描述符来估计对应关系，然后使用近似解来解决对应关系评估问题，将对应关系评估问题转换为gmm的概率密度。规定特征点xn和目标特征点ym之间的mlf描述符被认为是第n个高斯组件的质心与第m个数据的举例，因此获得的gmm概率密度函数为：

表示高斯混合模型的概率密度函数，pdf(ym|ρn)表示一个高斯组件ρn中对数据点ym的概率密度函数，其中包括n+1个高斯组件，m表示目标特征点集中的特征点，n表示源目标点集中的特征点，m表示目标特征点集中的特征点总数，n表示源目标点集中的特征点总数，每个高斯组件ρn表示源特征点集x中的xn，p(n)＝1/a，额外的第n+1个高斯组件用于消除冗余点的影响，使冗余点受参数ω约束(0＜ω＜1)；

一个高斯组件ρn中对数据点ym的概率密度如下：

pdf(ym|ρn)表示一个高斯组件ρn中对数据点ym的概率密度，mlf(t)表示尺度差异，σ²表示第一优化参数。

在获得由mlf引导的gmm的pdf之后，通过根据贝叶斯规则计算gmm的后验概率来估计的对应关系，其中gmm参数取前一次迭代的值：对应关系矩阵表达式如下：

最终获取由mlf描述符的一对多模糊对应关系p^old，记作矩阵p0，大小是n×m。同时获取相应的目标点集：

表示变换后的目标特征点集，y表示变换前的目标特征点集。在获得两个特征点集之间的对应关系矩阵之后，基于对应关系矩阵更新源特征点集x的空间变换。

由于源特征点集被构建为gmm，因此通过获得gmm的参数θ和σ²的最优解来完成空间变换更新。非刚性配准的空间变换更新可以被看做为一个参数优化过程，在本实施例中，通过最小化概率密度的负对数似然函数获取gmm的最优参数θ和σ²，进而获取最优参数矩阵。

概率密度的负对数似然函数表达式为：

本发明采用期望最大化(em)算法来解决这一优化问题。em算法包括计算后验概率(e步)和最大化期望(m步)这两步。根据之前获取的后验概率，将所述概率密度的负对数似然函数作为能量函数，展开如下：

通过最小化概率密度的负对数似然函数获取gmm的最优参数θ和σ²，σ²表示第一优化参数，θ表示第二优化参数，完成最大化期望。为了可以推导出能量函数的偏导数，本发明将其转换为纯矩阵形式。在此之前，首先基于高斯径向基函数(grbf)定义源特征点集x的空间变换c：

其中，w是高斯内核的n×d权重矩阵，g是由高斯径向基函数获得的n×n高斯内核矩阵，该矩阵中的某一元素的表达式为：

其中，n、m表示该元素在矩阵中的位置。w是高斯内核的n×d权重矩阵，优化参数θ即为w。能量函数的展开式如下：

用以下矩阵形式表示能量函数以进行偏导计算：

其中，trace(·)表示矩阵的迹，n×n矩阵px和m×m矩阵py是对角矩阵，分别由列向量p01和组成。1是一个列向量，其元素都是1。w是高斯内核的n×d权重矩阵，g是由高斯径向基函数获得的n×n高斯内核矩阵，x表示源特征点集，y表示目标特征点集。

通过求解以下偏导数的函数组来获得最佳参数：

关于全局几何约束，通过在能量函数中添加全局约束项以维持特征点集在空间变换更新时的全局结构稳定性，即采用一个正则化运算符以更新于全局结构稳定的非刚性变换。本实施例将此正则化运算符写作以下形式：

其中，λ是控制全局约束强度的权重参数。在能量函数的展开式添加全局约束项，得到第一能量函数，在第一能量函数中添加局部约束项huoqu第二能量函数。第一能量函数的表达式为：

qg(w,σ²)表示第一能量函数，q(w,σ²)表示能量函数，常量λ是控制全局约束强度的权重参数。

之后，得到新的偏导数组如下：

关于局部集合结构约束，本实施例基于s3中定义的局部特征描述符lgsf定义了局部几何结构约束项如下：

其中，η是控制局部约束强度的权重参数。约束的核心思想是通过在空间变换之前和之后判断x的局部相似性来约束空间变换中源特征集x的局部变形。权重参数η的值在本地约束项中具有重要影响。

然而，本实施例所提出的方法对局部约束项强度的要求因不同的情况和不同的配准过程而不同。因此，权重参数η的价值策略对所提方法的性能有一定影响。为了最大化方法的性能，在大多数情况下设计了对具有不同配准过程的局部约束项的约束强度的需求：在配准开始时，局部结构约束应用于空间变换。根据权重参数η的最大值设置源。配准过程中局部结构约束的约束强度逐渐变弱。因此，确定性退火技术用于定义权重参数η，如下所示：

其中，m是预先设置的权重参数η的最大值，t是当前的迭代次数，而常量c用于控制约束强度变化的速度。将局部几何结构约束项加到控制全局几何结构稳定性的能量函数中，使得该方法在更新空间变换的同时可以保持全局和局部稳定性。本文给出了局部结构约束的矩阵形式：

l＝ηtrace(d^td)；d＝uyy-ux(x+gw)

其中，i表示单位矩阵，y表示空间变换后的目标特征点集，表示空间变换前的目标特征点集。

假设特征点集t包括j个点，点tik是第k个特征点ti的最近相邻点，则此处的h(t)的定义为一个j×j的矩阵，该矩阵的每个元素定义如下：

即h(t)为一个稀疏矩阵，由ti的k个最近相邻点的权重组成。此定义通过设置配准前后点集的最近相邻点对应相同的点以加强约束项对点集局部结构稳定性的维护。将局部几何结构约束项加到第一能量函数中得到第二能量函数，第二能量函数表达式如下：

qgl(w,σ²)＝qg(w,σ²)+ηtrace(d^td)

由此得到新的公式：

由于每一项的推导：

因此，参数最优化的解为：

w＝(gpxg+2ησ²g^tu^tug)^-1(gp0y-gp0y)

px为n×n矩阵，py为m×m矩阵，px和py是对角矩阵，分别由列向量p01和组成，1是一个元素都是1的列向量。

s5、图像配准：基于后向方法的图像配准，根据在s4计算的源特征点集和目标特征点集之间的映射关系，获取变换后的图像。

具体地，在根据映射关系获得空间变换后的源特征点集后，即可与变换前的源点集x组合为然后实现图像配准。本发明采用反推法完成图像变换，即使用目标图像和图像上特征点作为控制点，根据点集的空间变换反推出变换后的图像。

通过将变换中的图像看做薄板，采用薄板样条(tps)变换模型完成空间变换，tps模型如下：

求得的是大小为(n+3)×3的矩阵，其中，是n×3的矩阵且它的第n行是ο是3×3的零矩阵，tps核是大小为n×n的矩阵。

根据矩形图像的逐像素索引获得一个矩形的网格点集的坐标其中网点集的数量为图像全部像素的数量z＝x(i^t)×y(i^t)。将网格δ^t作为tps模型的输入，得到具体如下：

和为两个矩阵。然后取其前两列作为变换后的网格点集的坐标其中大小为z×n的内核大小为z×3的矩阵的第z行为此时通过tps模型得到的网格点集即代表变换后的图像的位置，最后进行重采样获得变换后图像的内容。首先：

然后根据的坐标从待配准图像采样像素填充至变换后的图像除这些坐标外的位置的像素设置为0。在上述重采样过程中为了增强变换后的图像的平滑度，本发明在采样填充时使用了双三次插值法。

s6：获得变换后的图像，将其与参考图像进行像素级的融合得到融合图像，即完成多模态的视网膜图像融合。

首先，使用引导图像滤波来增强不同类型的视网膜图像的边缘。其次，采用多特征描述符来组合视网膜图像的边缘图和几何结构特征，以改进特征集的描述和排除冗余点。然后，多特征引导模型为特征集配准提供准确的引导，并且基于特征点的图像配准最终提供准确的图像配准。最后，基于准确的图像配准进行多模态的视网膜眼底图像融合。视网膜图像对的变化过程如说明书附图3、4、5所示，对于每组结果，第一行是源图像和目标图像，第二行是对应的边缘图，第三行是特征匹配结果，最后一行给出12×12的棋盘，用于交替显示转换后的图像和目标图像作为右侧的棋盘图像，变换后的图像显示在左侧。完整的流程示意图如说明书附图6所示。

本实施例提供了一种基于图像配准的多模态视网膜图像融合方法及系统。在高精度图像配准的基础上，进行有效的多模态眼底图像融合，便于眼科医生观察病变区域的变化，进一步协助相关视网膜疾病的诊断和治疗。

实施例2

本实施例在实施例1的基础上，将实施例1提出的一种基于图像配准的多模态视网膜图像融合方法模块化，形成一种基于图像配准的多模态视网膜图像融合系统，各模块示意图如说明书附图7所示。

一种基于图像配准的多模态视网膜图像融合系统，包括依次连接的图像输入单元1、图像预处理单元2、特征组合和提取单元3、特征点集配准单元4、图像配准单元5和图像融合单元6。

图像输入单元1：用于获取包括源图像和目标图像的视网膜图像对。每对视网膜图像对包括源图像和目标图像。

图像预处理单元2：用于对视网膜图像对进行预处理，获取视网膜边缘图像对，从每对视网膜边缘图像对中提取特征点集，特征点集包括从预处理后的源图像中提取的源特征点集和从预处理后的目标图像中提取的目标特征点集。预处理包括对比度增强、噪声消除等，引导图像滤波来增强不同类型的视网膜图像的边缘。

特征组合和提取单元3：用于组合多个特征描述符构建多特征差异描述子，引导对特征点集的特征提取。采用多特征描述符来组合视网膜图像的边缘图和几何结构特征，以改进特征集的描述和排除冗余点。

特征点集配准单元4：用于通过多特征差异描述子评估源特征点集与目标特征点集之间的对应关系，获取对应关系矩阵，根据对应关系矩阵更新源特征点集的空间变换。多特征引导模型为特征集配准提供准确的引导。

图像配准单元5：用于根据空间变换前后的源特征点集实现图像配准，结合目标图像获取空间变换后的图像。基于特征点的图像配准最终提供准确的图像配准。

图像融合单元6：用于将空间变换后的图像和预设参考图像进行像素级的融合得到融合图像。基于准确的图像配准进行多模态的视网膜眼底图像融合。

其中，特征点集配准单元4具体包括模型构建单元41、参数调整单元42和迭代单元43。具体包括：

模型构建单元41：用于在源特征点集上构造高斯混合模型，通过多特征差异描述子计算源特征点集与目标特征点集之间的对应关系矩阵。

参数调整单元42：用于基于对应关系矩阵，调整高斯混合模型的参数，基于全局几何结构约束和局部几何结构约束更新源特征点集的空间变化。

迭代单元43：用于迭代模型构建单元和参数调整单元，使源特征点集逐步接近目标特征点集，直至满足预设匹配关系。

其中多个特征描述符包括基于类sift尺度空间的边缘定向直方图描述符(eoh-sift描述符)、局部几何结构特征描述符、全局几何结构特征描述符。特征组合和提取单元3具体包括尺度差异单元31、局部特征单元32、全局特征单元33和差异构建单元34。具体包括：

尺度差异单元31：用于通过基于类sift尺度空间的边缘定向直方图，匹配位于相同场景、不同光谱条件下的视网膜边缘图像对上的特征点集，获取源特征点集和目标特征点集的尺度差异；

局部特征单元32：用于通过局部几何结构特征描述符，获取源特征点集和目标特征点集的局部特征差异矩阵；

全局特征单元33：用于通过全局几何结构特征描述符，获取源特征点集和目标特征点集的全局特征差异矩阵；

差异构建单元34：用于根据尺度差异、局部特征差异矩阵和全局特征差异矩阵构建多特征差异描述子。

本实施例在实施例1的基础上，提出了一种基于深度学习的样本图像生成系统，将实施例1的方法模块化，形成一种具体的系统，使其更具备实用性。

本发明针对现有技术，提供了一种基于图像配准的多模态视网膜图像融合方法与系统，解决了现有技术中的视网膜图像配准时匹配精度不足以及图像空间变换不准确的问题，采用本发明提供的方法能够获得多模态视网膜图像融合的最佳配准性能，并且在大多数情况下都优于目前最先进的方法。在高精度图像配准的基础上，进行有效的多模态眼底图像融合，便于眼科医生观察病变区域的变化，进一步协助相关视网膜疾病的诊断和治疗。

本领域普通技术人员应该明白，上述的本发明的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个计算装置上，或者分布在多个计算装置所组成的网络上，可选地，他们可以用计算机装置可执行的程序代码来实现，从而可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明不限制于任何特定的硬件和软件的结合。

注意，上述仅为本发明的较佳实施例及所运用技术原理。本领域技术人员会理解，本发明不限于这里的特定实施例，对本领域技术人员来说能够进行各种明显的变化、重新调整和替代而不会脱离本发明的保护范围。因此，虽然通过以上实施例对本发明进行了较为详细的说明，但是本发明不仅仅限于以上实施例，在不脱离本发明构思的情况下，还可以包括更多其他等效实施例，而本发明的范围由所附的权利要求范围决定。

以上公开的仅为本发明的几个具体实施场景，但是，本发明并非局限于此，任何本领域的技术人员能思之的变化都应落入本发明的保护范围。