CN109559292A - 基于卷积稀疏表示的多模态图像融合方法 - Google Patents

Abstract

本发明涉及一种基于卷积稀疏表示的多模态图像融合方法，首先，利用稀疏优化函数将源图像进行两尺度分解得到高频分量和低频分量；然后，将两尺度分解得到的高频和低频分量，根据多模态图像特点采用不同的融合策略，高频分量利用卷积稀疏表示对稀疏系数取最小值的融合策略，低频分量利用取平均的融合策略得到融合后图像的低频分量；最后将得到的融合后图像的高频分量和低频分量相加得到融合图像。相对其他三种融合方法，不论在主观视觉和客观评价指标上还是在计算效率上，本发明方法可以更好保留源图像的细节等纹理信息。

Description

基于卷积稀疏表示的多模态图像融合方法

技术领域

本发明涉及一种基于卷积稀疏表示的多模态图像融合方法，可以应用于各种民用的医学图像处理系统。

背景技术

随着计算机技术的高速发展和信息时代的到来，多样化的医学成像设备层出不穷，其应用贯穿于整个临床工作，为医学诊断和治疗提供多种模态医学图像，从解剖上和功能上多样化地描述人体形态信息或代谢信息。然而，受成像原理的限制，单一模态的图像都有各自的优缺点，很难全面地表达所有实用的信息，造成了诊断的局限性。若考虑对多幅非同源医学图像进行融合，整合和突出各自的互补信息，或能达到增强图像质量，减小冗余，为病灶的精准定位提供可靠的诊断治疗依据。因此多模态图像融合技术的研究有着深远的意义。近年，针对多模态图像融合技术，已有大量的图像融合算法被相继提出，并且迅速应用于多模态图像处理的各个方面，取得很好的效果。许多学者们已经提出多种方法进行解决。在多模态图像融合研究中，文献1“Burt P J,Adelson E H.The Laplacian Pyramid asa Compact Image Code[J].Readings in Computer Vision,1983,31(4):532-540.”和文献2“Li H,Manjunath B S,Mitra S K.Multisensor Image Fusion Using the WaveletTransform[C]Image Processing,1995,57(3):235-245.”提出了基于交叉双边滤波器、广义随机游动和马尔可夫随机场的多模态图像融合优化方法，虽然这些方法利用加权平均融合源图像来估计空域平滑和边缘对齐的权值，但是对于这种权值构造和全局优化方法可能会导致对比度降低、计算效率低下和权值的过平滑。此后为了减少建模负担和计算成本，学者们提出基于稀疏表示图像融合的方法，文献3“hang Q,Guo B L.Multifocus imagefusion using the nonsubsampled contourlet transform[J].Signal Processing,2009,89(7):1334-1346.”然而上述基于稀疏表示的算法一直是基于局部图像块的而不是基于整个图像块的，上述方法出现问题：图像细节保持能力有限，配准失调有很高的敏感性。

发明内容

要解决的技术问题

针对多模态图像融合后图像细节保持能力有限，配准失调有很高的敏感性等问题，本发明提出一种基于卷积稀疏表示的多模态图像融合方法。

技术方案

一种基于卷积稀疏表示的多模态图像融合方法，其特征在于步骤如下：

步骤1：多模态图像稀疏优化函数两尺度图像分解

首先将源图像I_n利用稀疏优化函数分解为一个低频分量和一个高频分量低频分量表示图像光谱信息，高频分量表示图像细节信息；低频部分利用优化数学表达式求解：

式中：表示源图像低频分量傅立叶变化，*表示卷积，η表示正则化参数，设置η为5，G_X和G_Y分别表示沿图像行和列的梯度滤波器；通过求解可得：

这样，通过求解逆傅立叶变化可以分别得到稀疏优化两尺度分解的低频分量和高频分量

步骤2：对于图像的高频部分其稀疏系数为c_n,m，m∈{1,2,...,M}，c_n,m通过求解如下卷积稀疏表示获得：

式中：D_m，G_l表示线性算子，表示一组M个字典滤波器集，{c_n,m}是一组稀疏系数集，{α_m}和{β_m}分别表示一组l₁和l₂范数的系数权重，是输入源图像的细节层，λ和μ是正则化参数；式中矩阵

那么式(3)的傅里叶变换表达式如下：

式中：D＝(D₀ D₁···)，β₀和β₁表示l₂范数的系数权重，α₀和α₁表示l₁范数的系数权重，分别为块矩阵；

引入一个辅助变量y₀，y₁，y₂，将式(4)变换为：

利用对偶变量，引入拉格朗日乘子u₀，u₁和u₂则通过迭代将式(5)的约束优化问题变为非约束优化问题形式：

u₀ ^(j+1)＝u₀ ^(j)+C^(j+1)-y₀ ^(j+1) (8)

u₁ ^(j+1)＝u₁ ^(j)+C^(j+1)-y₁ ^(j+1) (9)

u₂ ^(j+1)＝u₂ ^(j)+C^(j+1)-y₂ ^(j+1) (10)

式(7)通过下式给出：

S_γ(u)＝sign(u)⊙max(0,|u|-γ) (12)

由于式(6)的计算效率低下，在DFT域中用 和分别表示D_m、C_m、Γ₀、Γ₁、y₀、y₁、y₂、u₀、u₁、u₂和I^D；在DFT域中通过卷积定理，式(6)可变换为：

式中：

通过对式(13)进行二次优化的求解，对求偏导数并令其为0，可解得：

式中：矩阵为由M个N×N对角矩阵组成，M为滤波器的维数，N为源图像的维数，为一个MN×MN维的对称矩阵，和分别为对角矩阵；考虑到计算效率以及和ρI分别为对角矩阵，利用谢尔曼-莫里森公式对于式(14)求解，则式(14)利用谢尔曼-莫里森公式可求得：

然后对进行逆傅里叶变换，求得c_n,m；

假如c_n,1:M(x,y)表示c_n,m在空间域中位置(x,y)处的内容，c_n,1:M(x,y)是一个M维向量；根据基于稀疏表示的融合方法，利用c_n,1:M(x,y)的l₁范数作为源图像的活动水平度量；则活动水平图A_n(x,y)通过如下表达式获得：

A_n(x,y)＝||c_n,1:M(x,y)||₁ (16)

对A_n(x,y)利用基于窗口的平均策略以获得最终的活动水平图，可得：

式中：r决定窗口的大小，融合后的系数图为：

最后，融合图像的高频分量通过如下得到：

步骤3：多模态图像低频分量通过如下得到：

步骤4：将得到的融合后的高频分量和低频分量进行重构得到融合后的图像，则

有益效果

本发明提出一种卷积稀疏化表示的多模态图像融合算法，解决传统多模态融合算法在各自存在的缺陷。根据多模态图像融合的实验结果证明基于卷积稀疏化表示的多模态算法优于传统多模态图像融合算法的优点。

附图说明

图1本发明方法的基本流程图

图2多模态医学图像：(1)第一组源图像；(2)第二组源图像；(3)第三组源图像；(4)第四组源图像

图3第一组多模态源图像融合结果：(1)源图像1；(2)源图像2；(3)文献1；(4)文献2；(5)文献3；(6)本发明方法

图4第二组多模态源图像融合结果：(1)源图像1；(2)源图像2；(3)文献1；(4)文献2；(5)文献3；(6)本发明方法

图5第三组多模态源图像融合结果：(1)源图像1；(2)源图像2；(3)文献1；(4)文献2；(5)文献3；(6)本发明方法

图6第四组多模态源图像融合结果：(1)源图像1；(2)源图像2；(3)文献1；(4)文献2；(5)文献3；(6)本发明方法

具体实施方式

现结合实施例、附图对本发明作进一步描述：

用于实施的硬件环境是：实验环境为CPU Intel Core i5-5200U@2.20GHz，内存为4GB，采用MATLAB R2014a编程。本发明采用多模态图像进行融合处理。

本发明方法的基本流程如图1所示，实验源图像数据如图2所示，具体实施如下：

步骤一：多模态图像稀疏优化函数两尺度图像分解

首先将源图像I_n利用稀疏优化函数分解为一个低频分量和一个高频分量分解的图像的低频部分主要表征图像光谱信息；图像的高频部分主要表示图像细节信息。其中低频部分利用优化数学表达式求解：

式中：表示源图像低频分量傅立叶变化，*表示卷积，η表示正则化参数，设置η为5，G_X和G_Y分别表示沿图像行和列的梯度滤波器。通过求解可得：

这样，通过求解逆傅立叶变化可以得到稀疏优化两尺度分解的低频分量高频分量可通过源图像I_n减去低频分量获得。

步骤二：多模态图像的高频分量

多模态图像的高频分量反映图像边缘、纹理等细节信息。对于图像的高频部分其稀疏系数为c_n,m，m∈{1,2,...,M}，c_n,m可以通过求解如下卷积稀疏表示获得：

式中：D_m，G_l表示线性算子，表示一组M个字典滤波器集，{c_n,m}是一组稀疏系数集，{α_m}和{β_m}分别表示一组l₁和l₂范数的系数权重，是输入源图像的细节层，λ和μ是正则化参数。式中矩阵

式中：G_l表示线性算子，β₀和β₁表示l₂范数的系数权重，那么式(3)的傅里叶变换表达式如下：

式中：D＝(D₀ D₁···)，α₀和α₁表示l₁范数的系数权重，分别为块矩阵。引入一个辅助变量y₀，y₁，y₂，将式(4)变换为：

利用对偶变量，引入拉格朗日乘子u₀，u₁和u₂则通过迭代将式(5)的约束优化问题变为非约束优化问题形式：

u₀ ^(j+1)＝u₀ ^(j)+C^(j+1)-y₀ ^(j+1) (8)

u₁ ^(j+1)＝u₁ ^(j)+C^(j+1)-y₁ ^(j+1) (9)

u₂ ^(j+1)＝u₂ ^(j)+C^(j+1)-y₂ ^(j+1) (10)

式(7)通过下式给出：

S_γ(u)＝sign(u)⊙max(0,|u|-γ) (12)

由于式(6)的计算效率低下，在DFT域中用 和分别表示D_m、C_m、Γ₀、Γ₁、y₀、y₁、y₂、u₀、u₁、u₂和I^D。在DFT域中通过卷积定理，式(6)可变换为：

式中：

通过对式(13)进行二次优化的求解，对求偏导数并令其为0，可解得：

式中：矩阵为由M个N×N对角矩阵组成，M为滤波器的维数，N为源图像的维数，为一个MN×MN维的对称矩阵，和分别为对角矩阵。考虑到计算效率以及和ρI分别为对角矩阵，利用谢尔曼-莫里森公式对于式(14)求解，则式(14)利用谢尔曼-莫里森公式可求得：

然后对进行逆傅里叶变换，求得c_n,m。

假如c_n,1:M(x,y)表示c_n,m在空间域中位置(x,y)处的内容，c_n,1:M(x,y)是一个M维向量。根据基于稀疏表示的融合方法，利用c_n,1:M(x,y)的l₁范数作为源图像的活动水平度量。则活动水平图A_n(x,y)通过如下表达式获得：

A_n(x,y)＝||c_n,1:M(x,y)||₁ (16)

对A_n(x,y)利用基于窗口的平均策略以获得最终的活动水平图，可得：

式中：r决定窗口的大小，在多模态图像融合中，由于源图像的小尺度细节往往存在，本发明采用“较小的r”融合策略完成多模态图像高频分量融合，融合后的系数图为：

最后，融合图像的高频分量通过如下得到：

步骤三：多模态图像的低频分量

根据多模态图像特点，本发明采用“平均的融合策略”完成多模态低频分量融合，保留源图像的纹理等细节信息，满足人眼视觉特性，多模态图像低频分量通过如下得到：

步骤四：对两尺度图像的高频和低频分量进行重构，得到多模态融合图像

将得到的融合后的高频分量和低频分量进行重构得到融合后的图像，则

下面结合图1～图6、表1对本发明的效果做进一步描述。

图2是四组多模态源图像。

1.实验条件

实验环境为CPU Intel Core i5-5200U@2.20GHz，内存为4GB，采用MATLAB R2014a编程。本发明采用多模态图像集(256×256)如图2所示。

2.实验内容

图3-图6是多模态图像和融合后图像的对比图。

将四组CT/MRI多模态医学源图像分别经过基于文献1、文献2、文献3及本发明方法进行融合处理，分别对各算法的融合图像进行实验分析。图3-图6为四组CT/MRI多模态医学源图像及其不同算法融合结果图。

从主观视觉上分析，对于融合图像，文献1方法虽然和别的方法相比，在亮度和对比度上优于别的方法，但是局部细节信息丢失严重，(如图3(3)～图6(3))；文献2方法出现严重的失真和细节信息丢失现象，(如图3(4)～图6(4))；文献3方法所得的融合图像不仅丢失大量的图像细节信息，而且还出现严重的伪影现象，(如图3(5)～图6(5))；通过以上几种融合算法的效果对比，本发明方法融合图像主观效果优于其他算法，所得的融合图像不仅可以保留源图像的大量细节信息，又不会产生可见的伪影和亮度失真。

从客观定量分析，采用图像互信息MI、信息结构相似度Q_Y、峰值信噪比PSNR和边缘保持度Q^AB/F等客观评价指标，对以上几种融合方法效果进行定量评价。对于四组医学图像，本发明融合的方法与其他算法的定量评价比较见表1。从表1可以看到在第一组、第三组和第四组中本发明融合方法的MI、Q^AB/F、Q_Y和SD都是最大；在第二组中本发明融合方法的Q^AB/F、SD和MI，只有文献3的Q_Y值最大，但Q_Y值只是略微相差一点。

表1多模态医学图像不同融合方法的客观评价指标

总之，在这四组客观数据比较中，本发明方法的指标基本上都优于文献的方法，综合主观视觉和客观定量评价指标，本发明方法更优于其它算法，可以更有效地从多模态源图像中提取出细节信息，并凸显出显著性特征。

Claims (1)

1.一种基于卷积稀疏表示的多模态图像融合方法，其特征在于步骤如下：

步骤1：多模态图像稀疏优化函数两尺度图像分解

首先将源图像I_n利用稀疏优化函数分解为一个低频分量和一个高频分量低频分量表示图像光谱信息，高频分量表示图像细节信息；低频部分利用优化数学表达式求解：

式中：表示源图像低频分量傅立叶变化，*表示卷积，η表示正则化参数，设置η为5，G_X和G_Y分别表示沿图像行和列的梯度滤波器；通过求解可得：

这样，通过求解逆傅立叶变化可以分别得到稀疏优化两尺度分解的低频分量和高频分量

步骤2：对于图像的高频部分其稀疏系数为c_n,m，m∈{1,2,...,M}，c_n,m通过求解如下卷积稀疏表示获得：

式中：D_m，G_l表示线性算子，表示一组M个字典滤波器集，{c_n,m}是一组稀疏系数集，{α_m}和{β_m}分别表示一组l₁和l₂范数的系数权重，是输入源图像的细节层，λ和μ是正则化参数；式中矩阵

那么式(3)的傅里叶变换表达式如下：

式中：D＝(D₀ D₁…)，β₀和β₁表示l₂范数的系数权重，α₀和α₁表示l₁范数的系数权重，分别为块矩阵；

引入一个辅助变量y₀，y₁，y₂，将式(4)变换为：

利用对偶变量，引入拉格朗日乘子u₀，u₁和u₂则通过迭代将式(5)的约束优化问题变为非约束优化问题形式：

u₀ ^(j+1)＝u₀ ^(j)+C^(j+1)-y₀ ^(j+1) (8)

u₁ ^(j+1)＝u₁ ^(j)+C^(j+1)-y₁ ^(j+1) (9)

u₂ ^(j+1)＝u₂ ^(j)+C^(j+1)-y₂ ^(j+1)(10)

式(7)通过下式给出：

S_γ(u)＝sign(u)⊙max(0,|u|-γ) (12)

由于式(6)的计算效率低下，在DFT域中用 和分别表示D_m、C_m、Γ₀、Γ₁、y₀、y₁、y₂、u₀、u₁、u₂和I^D；在DFT域中通过卷积定理，式(6)可变换为：

式中：

通过对式(13)进行二次优化的求解，对求偏导数并令其为0，可解得：

式中：矩阵为由M个N×N对角矩阵组成，M为滤波器的维数，N为源图像的维数，为一个MN×MN维的对称矩阵，和分别为对角矩阵；考虑到计算效率以及和ρ^I分别为对角矩阵，利用谢尔曼-莫里森公式对于式(14)求解，则式(14)利用谢尔曼-莫里森公式可求得：

然后对进行逆傅里叶变换，求得c_n,m；

假如c_n,1:M(x,y)表示c_n,m在空间域中位置(x,y)处的内容，c_n,1:M(x,y)是一个M维向量；根据基于稀疏表示的融合方法，利用c_n,1:M(x,y)的l₁范数作为源图像的活动水平度量；则活动水平图A_n(x,y)通过如下表达式获得：

A_n(x,y)＝||c_n,1:M(x,y)||₁ (16)

对A_n(x,y)利用基于窗口的平均策略以获得最终的活动水平图，可得：

式中：r决定窗口的大小，融合后的系数图为：

最后，融合图像的高频分量通过如下得到：

步骤3：多模态图像低频分量通过如下得到：

步骤4：将得到的融合后的高频分量和低频分量进行重构得到融合后的图像，则