CN108596866B - 一种基于稀疏低秩分解和视觉显著性结合的医学图像融合方法 - Google Patents

Abstract

本发明涉及一种基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，属于数字图像处理技术领域。本发明设计了一种有效的图像低秩稀疏成分分解模型。其中，分解得到的低秩成分反映了图像的整体轮廓和亮度信息，稀疏成分反映了图像的高频细节信息。本明针对图像不同的成分融合，提出了不同的融合方案。在融合低秩成分时，采用传统的“绝对值”取大的融合策略，以保留源图像更多的亮度信息；而对于稀疏成分，本申请提出一种视觉显著性度量的方法来保留更多的显著性稀疏成分。本发明提出的图像融合的方法无论是从人眼视觉上还是从客观评价指标上，本文方法均优于其它几种融合方法。

Description

一种基于稀疏低秩分解和视觉显著性结合的医学图像融合 方法

技术领域

本发明涉及一种基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，属于数字图像处理技术领域。

背景技术

医学图像融合是医学图像处理技术的重要分支，其目的是把来自不同成像模式的多幅源图像根据某种方式进行综合形成一幅图像，从而实现对目标物的精确、全面描述。目前，多模态医学图像融合技术已被成功应用于医疗诊断中。

为获得视觉效果更好的融合图像，研究者提出了多种医学图像融合方法。在这些方法中，最为经典是多尺度变换域的融合方法。这类方法通常先通过某种变换，将输入的图像进行分解，对分解后的结果进行融合处理，最后再经过逆变换来获得最终的融合图像。常见的变换方法有金字塔变换、离散小波变换，非采样变换等。然而，基于小波变换的融合方法，由于其不能实现对图像的稀疏表示，从而影响最终融合结果的视觉效果；同时，大部分小波变换不具有平移不变性，融合结果中容易引入一些虚假信息。

相对于上述多尺度变换方法，基于稀疏表示的医学图像融合能有效克服多尺度变换方法带来的一些不足，因此近几年受到了广泛关注。在基于稀疏表示的图像融合中，字典的构建是影响融合结果的关键因素之一。通常，字典可通过解析的方式来产生，这类字典包含常见的DCT字典、wavelet字典、Curvelet字典等。然而，这种基于解析方式设计的字典很难实现对自然图像复杂结构的稀疏刻画，也很难对图像内容进行自适应地描述。相反，通过某种算法如KSVD从训练样本中学习得到的字典可有效避免解析式构成的字典所具有的缺陷，能实现对信号或图像结构的自适应描述。

为得到表达能力更强的字典，进一步提升融合图像的视觉效果，研究者通过字典的设计实现了融合结果视觉效果的提升，但是对图像的不同成分，使用的是同一个字典来对图像内容进行表示。事实上，图像的内容可以看成具有不同空间形态结构成分构成的。采用固定统一的字典很难实现对图像不同成分的有效表示。虽然，已有文献注意到此问题，并提出采用不同字典来对图像不同的成分进行稀疏表示，但其所使用的字典是通过解析方式设计得到的，很难实现对图像内容的自适应描述。同时，现有方法并未考虑到图像不同成分的分解和表述不仅与字典有关，更与图像分解所采用的模型有关。此外，传统基于稀疏表示的医学图像融合方法，大都采用“绝对值取大”的融合策略来获取融合图像的稀疏表示系数，从而忽略了图像的对比度信息。

发明内容

本发明提供了一种基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，以用于保留源图像的亮度信息及保留显著性特征的基础上获得较好的融合图像。

本发明的技术方案是：一种基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，所述方法的具体步骤如下：

Step1、收集多幅多源图像，以t*t为大小的滑动窗口从多源图像中采集N个图像块作为训练样本的数据，采集到的样本数据记为Y；

Step2、通过字典学习模型对采集到的样本数据Y进行训练，得到低秩字典D_l和稀疏字典D_s；

Step3、从Step2中得到低秩字典D_l、稀疏字典D_s后，取待融合图像A、B两张，以t*t为大小的滑动窗口从待融合图像A、B中分别采集N个图像块数据，每个图像块数作为一个列向量，再将列向量组成矩阵，图A、图B采集到的数据组成的矩阵分别为X_A、X_B；然后用鲁棒主成分分析法将X_A、X_B分解为低秩、稀疏成分；再利用字典D_l、D_s通过低秩稀疏分解模型求解得到待融合图像A和待融合图像B对应的低秩、稀疏成分的编码系数β_l、β_s；

Step4、根据Step3中得到的待融合图像A和待融合图像B对应的低秩成分的编码系数β_l，采用编码系数的l₁范数值最大的方法，得到融合图像的低秩成分的编码系数

再由待融合图像A和待融合图像B对应的稀疏成分的编码系数β_s，引入图像对比度，并结合编码系数的l₁范数值最大的方法计算得到融合图像的稀疏成分的编码系数

Step5、通过Step4计算得到的融合图像的低秩、稀疏成分的编码系数

和

以及Step2中的得到的低秩、稀疏成分的字典D_l,D_s可得到融合图像的低秩、稀疏成分

将两部分相加，可得到融合图像

所述步骤Step2中，字典学习模型如下：

式中，Y＝[y₁,y₂,…,y_N]∈R^M，y_i∈R^M(i＝1,2,…,N)为Y中第i列，为滑动窗口采集的训练样本图像的数据，R^M为向量空间，M是向量空间的维数，N是图像块的个数，Y_l、Y_s分别为用RPCA分解得到的低秩、稀疏成分矩阵，D_l∈R^M×K、D_s∈R^M×K分别代表M行K列的低秩字典和稀疏字典，M〈〈K，α_l＝[α_l,1,α_l,2,…,α_l,N]∈R^K×N，α_s＝[α_s,1,α_s,2,…,α_s,N]∈R^K×N分别为低秩、稀疏分量对应的编码系数，α_l,i表示低秩分量在低秩字典D_l下的第i个原子对应的稀疏编码向量，α_s,i表示稀疏分量在稀疏字典D_s下的第i个原子对应的稀疏编码向量，μ₁、μ₂、λ₁、λ₂为平衡参数，||·||_F表示F范数运算符，||·||_*表示核范数运算符，||·||₁表示l₁范数，

表示F范数平方运算符。

所述步骤Step3中，低秩稀疏分解模型如下所示：

式中，X为滑动窗口采集的待融合图像的数据作为列向量组成的矩阵，X＝[x₁,x₂,…,x_N]∈R^M，x_i(i＝1,2,…,N)为X中第i列，R^M为向量空间，M是向量空间的维数，N是图像块的个数；D_l∈R^M×K、D_s∈R^M×K分别为Step2中训练得到M行K列的低秩、稀疏字典，β_l＝[β_l,1,β_l,2,…,β_l,N]，β_s＝[β_s,1,β_s,2,…,β_s,N]，β_l∈R^K×N，β_s∈R^K×N分别表示待融合图像K行N列低秩、稀疏成分的编码系数，β_i,l表示低秩分量在字典D_l下的第i个原子对应的稀疏编码向量，β_i,s表示稀疏分量在字典D_s下的第i个原子对应的稀疏编码向量，μ₁、μ₂、λ₁、λ₂为平衡参数；T表示矩阵的转置，||·||_F表示F范数运算符，||·||_*表示核范数运算符，||·||₁表示l₁范数，

表示F范数平方运算符。

所述低秩稀疏分解模型求解如下：引入辅助变量X_l、X_s，使得X_l＝D_lβ_l、X_s＝D_sβ_s，则公式(2)可转化为公式(3)：

使用交替迭代算法，即更新X_l时固定X_s、β_l、β_s不变，目标函数如下所示：

同理更新X_s时固定X_l、β_l、β_s不变，目标函数如下所示：

在更新后的X_l和X_s的基础上，得到更新低秩、稀疏成分的编码系数β_l和β_s的目标函数为：

采用正交匹配追踪算法求解公式(6)、(7)，能得到待融合图像A、B的低秩、稀疏成分的编码系数。

所述步骤Step4中，编码系数的l₁范数值最大的方法如下：

公式(8)中，

表示融合图像的低秩成分的编码系数的第n＝1,2,…,N列，

表示第r＝1，2,...,P张待融合图像的低秩成分的编码系数β_l的第n＝1,2,…,N列，是一个列向量，则

表示的是融合图像的低秩成分的编码系数；N是图像块的个数。

所述步骤Step4中，融合图像的稀疏成分的编码系数计算如下：

公式(9)-(10)中，

表示融合图像的稀疏成分的编码系数的第n＝1,2,…,N列，

表示第j＝1,2,...,P,P+1张待融合图像的稀疏成分的编码系数β_s的第n＝1,2,…,N列，

表示第k＝1，2,...,P张待融合图像的稀疏成分的编码系数β_s的第n＝1,2,…,N列，是一个列向量，并且是稀疏成分的编码系数的列中l₁范数最大的一个；

表示的是融合图像的稀疏成分的编码系数，

表示第j＝1,2,...,P,P+1张待融合图像的第n列的低秩、稀疏成分对比度的l₁范数最大的一个；N是图像块的个数。

本发明的有益效果是：

1、本发明设计了一种有效的图像低秩稀疏成分分解模型。其中，分解得到的低秩成分反映了图像的整体轮廓和亮度信息，稀疏成分反映了图像的高频细节信息。

2、本发明针对图像不同的成分融合，提出了不同的融合方案。在融合低秩成分时，采用传统的“绝对值”取大的融合策略，以保留源图像更多的亮度信息；而对于稀疏成分，本文提出一种视觉显著性度量的方法来保留更多的显著性稀疏成分。

3、本发明提出的图像融合的方法无论是从人眼视觉上还是从客观评价指标上，本文方法均优于其它几种融合方法。

附图说明

图1为本发明的流程图；

图2为6张训练样本图像；

图3为本发明的CT和MTR待融合源图像；

图4为NSCT、SR、Kim、NSCT-SR、Zhu与本发明方法的融合图像对比。

具体实施方式

实施例1：如图1-4所示，一种基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，所述方法的具体步骤如下：

Step1、首先从互联网上收集6张细节信息丰富的多源图像，如图2所示，其中包括医学图像、人物图像和建筑物图像(多源图像为三幅以上，至少包括一幅医学图像、一幅人物图像和一幅建筑物图像)，以t*t为大小的滑动窗口从多源图像中采集N个图像块作为训练样本的数据，采集到的样本数据记为Y；

Step2、通过字典学习模型对采集到的样本数据Y进行训练，得到低秩字典D_l和稀疏字典D_s。

所述字典学习模型如下：

式(1)中，Y＝[y₁,y₂,…,y_N]∈R^M，Y为滑动窗口采集的训练样本图像的数据，R^M为向量空间，M是向量空间的维数，N是图像块的个数。其中y_i∈R^M(i＝1,2,…,N)为Y中第i列。Y_l、Y_s分别为用RPCA(鲁棒主成分分析法)分解得到的低秩、稀疏成分矩阵，D_l∈R^M×K，D_s∈R^M×K(其中M的值远小于K)，分别代表M行K列的低秩字典和稀疏字典，α_l＝[α_l,1,α_l,2,…,α_l,N]∈R^{K ×N}，α_s＝[α_s,1,α_s,2,…,α_s,N]∈R^K×N分别为低秩、稀疏分量对应的编码系数。α_l,i和α_s,i分别表示低秩分量和稀疏分量在字典D_l和D_s下的第i个原子对应的稀疏编码向量，μ₁、μ₂、λ₁、λ₂为平衡参数，其中λ₁＝λ₂＝0.01，μ₁＝μ₂＝1。||·||_F表示F范数运算符，||·||_*表示核范数运算符，||·||₁表示l₁范数。

表示F范数平方运算符。

Step3、从Step2中得到低秩D_l、稀疏字典D_s后，输入两幅像素大小为256×256的CT和MRI图像(如图3(a)、3(b))，对CT和MRI图像以滑动窗口的形式分别采集两张图像的数据，每个窗口采集到的数据作为一个列向量，再将列向量组成矩阵。CT和MRI图像采集到的矩阵分别为X_A，X_B，然后用RPCA算法分别将X_A,X_B分解为低秩、稀疏成分。再利用字典D_l,D_s通过低秩、稀疏分解模型求解得到图像A和图像B的对应的低秩、稀疏表示系数β_l,β_s。

所述步骤Step3中，低秩稀疏分解模型如下所示：

式中X为滑动窗口采集的待融合图像的数据作为列向量组成的矩阵，X＝[x₁,x₂,…,x_N]∈R^M，R^M为向量空间，M是向量空间的维数，N是图像块的个数。其中x_i(i＝1,2,…,N)为X中第i列。D_l∈R^M×K，D_s∈R^M×K分别为Step2中训练得到M行K列的低秩、稀疏字典，β_l＝[β_l,1,β_l,2,…,β_l,N]，β_s＝[β_s,,1,β_s,2,…,β_s,N]，β_l∈R^K×N，β_s∈R^K×N分别表示待融合图像K行N列低秩、稀疏成分的编码系数。β_i,l和β_i,s分别表示低秩分量和稀疏分量在字典D_l和D_s下的第i个原子对应的稀疏编码向量，λ₁＝λ₂＝0.01，μ₁＝μ₂＝1。其中T表示矩阵的转置，||·||_F表示F范数运算符，||·||_*表示核范数运算符，||·||₁表示l₁范数。

表示F范数平方运算符。

所述低秩稀疏分解模型求解如下：引入辅助变量X_l、X_s，使得X_l＝D_lβ_l、X_s＝D_sβ_s，则公式(2)可转化为公式(3)：

使用交替迭代算法，即更新X_l时固定X_s、β_l、β_s不变，目标函数如下所示：

同理更新X_s时固定X_l、β_l、β_s不变，目标函数如下所示：

在更新后的X_l和X_s的基础上，得到更新低秩、稀疏成分的编码系数β_l和β_s的目标函数为：

采用正交匹配追踪算法求解公式(6)、(7)，能得到图像A、B的低秩、稀疏成分的编码系数。通过上述方法即可分别求解得到CT图像的的低秩、稀疏成分的编码系数β_Al,β_As，MRI图像的低秩、稀疏成分的编码系数β_Bl,β_Bs。

Step4、根据step3中得到CT和MRI图像的低秩成分的编码系数β_Al,β_Bl，采用编码系数的l₁范数值最大的方法，也就是“绝对值”取大的融合策略，得到融合图像的低秩成分的编码系数

再由CT和MRI图像的稀疏成分的编码系数β_As,β_Bs，引入图像对比对度L的值，并结合“绝对值”最大策略计算得到融合图像的稀疏成分的编码系数

上述“绝对值”最大策略计算方式如下：

公式(8)中，

表示融合图像低秩成分的编码系数的第n＝1,2,…,N列，

表示第r＝1，2,...,P张待融合图像的低秩成分的编码系数β_l的第n＝1,2,…,N列，是一个列向量，则

表示的是融合图像的低秩成分的编码系数；N是图像块的个数。

所述步骤Step4中，融合图像的稀疏成分的编码系数计算如下：

公式(9)-(10)中，

表示融合图像的稀疏成分的编码系数的第n＝1,2,…,N列，

表示第j＝1,2,...,P,P+1张待融合图像的稀疏成分的编码系数β_s的第n＝1,2,…,N列，

表示第k＝1，2,...,P张待融合图像的稀疏成分的编码系数β_s的第n＝1,2,…,N列，是一个列向量，并且是的稀疏成分的编码系数的列中l₁范数最大的一个；

表示的是融合图像的稀疏成分的编码系数，

表示第j＝1,2,...,P,P+1张待融合图像的第n列的低秩、稀疏成分对比度的l₁范数最大的一个；N是图像块的个数。

Step5、通过Step4计算得到的融合图像的低秩、稀疏成分的编码系数

和

以及Step2中的得到的低秩、稀疏成分的字典D_l,D_s可得到融合图像的低秩、稀疏成分

将两部分相加，可得到融合图像X_F，如公式(11)所示。用matlab函数将数据矩阵X_F重新排成图像，可得到将CT和MRI的图像融合后的结果图，如图4(f)，即为本发明的方法的融合结果。

另外图4的(a)～(e)分别表示NSCT，SR，Kim，NSCT-SR，Zhu的方法融合得到的结果图。由图可知：基于NSCT的方法融合图像整体对比度较低。基于SR的方法对于图像的边缘信息细节信息的保留，但是不如本文方法和Zhu、Kim以及NSCT-SR的方法。Zhu以及Kim等提出的方法表现出了较好好的融合性能，这主要取决于两种方法均采用了信息量更为丰富的字典来进行稀疏表示。但和本发明方法图4(f)相比仍略显不足，这不仅仅是不同方法对不同成分采用了不同的字典进行稀疏表示，还得益于本文在分解图像时设计了更为合理的图像分解模型。

所述步骤Step1中，用滑动窗口大小为t×t(t＝8)，滑动窗口采集的数据形成的256维的矩阵就是训练字典的样本数据；

所述步骤Step2中，更新典以原子为单位更新，更新稀疏成分的编码系数以列为单位更新，直到所有列和原子更新完成，得到最后的低秩、稀疏字典。

所述步骤Step3中，由于低秩和稀疏成分具有不同的空间形态结构，本发明设计了低秩与稀疏成分分解模型，并通过低秩、稀疏字典，对不同成分进行了低秩、稀疏表达。

所述步骤Step4中，对图像不同的成分融合，提出了不同的融合方案。融合后图像的大小与待融合图像的大小要严格一致(均为256×256)，并把重构的图像以二值图像的方式显示。

除主观上对不同融合方法得到融合结果进行评价之外，客观评价也是必不可少的。本发明性能的客观评价指标包括Q_IM(MI,Mutual Information)、Q_TE(TE,TsallisEntropy)、Q_NCIE(NCIE,Nonlinear Correlation Information Entropy)、Q_S(PM,PiellaMetric),Q_G(GM,Gradient Based Metric)和Q_P(Metric-based on PhaseCongruency)。Q_IM和Q_TE可以用来度量源图像中有多少信息转移到了融合图像中；Q_NCIE可度量融合图像与源图像信息之间的相关性大小；Q_S用来衡量源图像块A相对于源图像块B的重要性；Q_G用来度量融合图像与源图像之间细节信息变化大小；Q_P度量融合图像中包含多少源图像的纹理和边缘信息。各个指标值越大，融合质量就更高。

表1给出了CT和RMI图像不同融合方法的客观评价数据，从数据中可以发现本发明方法具有最好的客观评价值，因此本文提出的方法相对于其它方法具有更优的融合性能。

表1客观质量评估不同方法的融合效果(图3.(a)-(f))

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

Claims (4)

1.一种基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，其特征在于：所述方法的具体步骤如下：

Step1、收集多幅多源图像，以t*t为大小的滑动窗口从多源图像中采集N个图像块作为训练样本的数据，采集到的样本数据记为Y；

Step2、通过字典学习模型对采集到的样本数据Y进行训练，得到低秩字典D_l和稀疏字典D_s；

Step3、从Step2中得到低秩字典D_l、稀疏字典D_s后，取待融合图像A、B两张，以t*t为大小的滑动窗口从待融合图像A、B中分别采集N个图像块数据，每个图像块数作为一个列向量，再将列向量组成矩阵，图A、图B采集到的数据组成的矩阵分别为X_A、X_B；然后用鲁棒主成分分析法将X_A、X_B分解为低秩、稀疏成分；再利用字典D_l、D_s通过低秩稀疏分解模型求解得到待融合图像A和待融合图像B对应的低秩、稀疏成分的编码系数β_l、β_s；

Step4、根据Step3中得到的待融合图像A和待融合图像B对应的低秩成分的编码系数β_l，采用编码系数的l₁范数值最大的方法，得到融合图像的低秩成分的编码系数

再由待融合图像A和待融合图像B对应的稀疏成分的编码系数β_s，引入图像对比度，并结合编码系数的l₁范数值最大的方法计算得到融合图像的稀疏成分的编码系数

Step5、通过Step4计算得到的融合图像的低秩、稀疏成分的编码系数

和

以及Step2中的得到的低秩、稀疏成分的字典D_l,D_s可得到融合图像的低秩、稀疏成分

将两部分相加，可得到融合图像

所述步骤Step2中，字典学习模型如下：

式中，Y＝[y₁,y₂,…,y_N]∈R^M，y_i∈R^M，i＝1,2,...,N为Y中第i列，为滑动窗口采集的训练样本图像的数据，R^M为向量空间，M是向量空间的维数，N是图像块的个数，Y_l、Y_s分别为用RPCA分解得到的低秩、稀疏成分矩阵，D_l∈R^M×K、D_s∈R^M×K分别代表M行K列的低秩字典和稀疏字典，M〈〈K，α_l＝[α_l,1,α_l,2,…,α_l,N]∈R^K×N，α_s＝[α_s,1,α_s,2,…,α_s,N]∈R^K×N分别为低秩、稀疏分量对应的编码系数，α_l,i表示低秩分量在低秩字典D_l下的第i个原子对应的稀疏编码向量，α_s,i表示稀疏分量在稀疏字典D_s下的第i个原子对应的稀疏编码向量，μ₁、μ₂、λ₁、λ₂为平衡参数，||·||_F表示F范数运算符，||·||_*表示核范数运算符，||·||₁表示l₁范数，

表示F范数平方运算符；

所述步骤Step3中，低秩稀疏分解模型如下所示：

式中，X为滑动窗口采集的待融合图像的数据作为列向量组成的矩阵，X＝[x₁,x₂,…,x_N]∈R^M，x_i∈R^M，i＝1,2,...,N为X中第i列，R^M为向量空间，M是向量空间的维数，N是图像块的个数；D_l∈R^M×K、D_s∈R^M×K分别为Step2中训练得到M行K列的低秩、稀疏字典，β_l＝[β_l,1,β_l,2,…,β_l,N]，β_s＝[β_s,1,β_s,2,…,β_s,N]，β_l∈R^K×N，β_s∈R^K×N分别表示待融合图像K行N列低秩、稀疏成分的编码系数，β_i,l表示低秩分量在字典D_l下的第i个原子对应的稀疏编码向量，β_i,s表示稀疏分量在字典D_s下的第i个原子对应的稀疏编码向量，μ₁、μ₂、λ₁、λ₂为平衡参数；T表示矩阵的转置，||·||_F表示F范数运算符，||·||_*表示核范数运算符，||·||₁表示l₁范数，

表示F范数平方运算符。

2.根据权利要求1所述的基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，其特征在于：所述低秩稀疏分解模型求解如下：引入辅助变量X_l、X_s，使得X_l＝D_lβ_l、X_s＝D_sβ_s，则公式(2)可转化为公式(3)：

使用交替迭代算法，即更新X_l时固定X_s、β_l、β_s不变，目标函数如下所示：

同理更新X_s时固定X_l、β_l、β_s不变，目标函数如下所示：

在更新后的X_l和X_s的基础上，得到更新低秩、稀疏成分的编码系数β_l和β_s的目标函数为：

采用正交匹配追踪算法求解公式(6)、(7)，能得到待融合图像A、B的低秩、稀疏成分的编码系数。

3.根据权利要求1所述的基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，其特征在于：所述步骤Step4中，编码系数的l₁范数值最大的方法如下：

公式(8)中，

表示融合图像的低秩成分的编码系数的第n＝1,2,…,N列，

表示第r＝1,2,...,P张待融合图像的低秩成分的编码系数β_l的第n＝1,2,…,N列，是一个列向量，则

表示的是融合图像的低秩成分的编码系数；N是图像块的个数。

4.根据权利要求1所述的基于稀疏低秩分解和视觉显著性结合的医学图像融合方法，其特征在于：所述步骤Step4中，融合图像的稀疏成分的编码系数计算如下：

公式(9)中，

表示融合图像的稀疏成分的编码系数的第n＝1,2,…,N列，

表示第j＝1,2,...,P张待融合图像的稀疏成分的编码系数β_s的第n＝1,2,…,N列；

表示的是融合图像的稀疏成分的编码系数，

表示第j＝1,2,...,P张待融合图像的第n列的低秩、稀疏成分对比度的l₁范数，

表示第j＝1,2,...,P张待融合图像的低秩成分的编码系数β_l的第n＝1,2,…,N列；N是图像块的个数。

Images