CN102968781B - 基于nsct和稀疏表示的图像融合方法 - Google Patents

Abstract

本发明提供了一种基于NSCT和稀疏表示的图像融合方法，对稀疏度较差的低频子带系数学***自适应调整权重融合；对稀疏度较高的高频方向子带系数采用同一尺度下方向子带绝对值和取最大的方法融合，以捕获源图像中的显著特征，最终提高融合效果。

Description

基于NSCT和稀疏表示的图像融合方法

技术领域

本发明涉及一种图像融合方法。

背景技术

近年来，基于非下采样Contourlet变换(Non-Subsampled Contourlet Transform，NSCT)以其具有平移不变、多分辨率、多方向和各向异性的图像表示能力，并且能有效克服传统小波变换不能处理2D或更高维奇异性的问题，成功用于图像融合领域并取得较优的融合效果。然而，在图像融合问题中，我们希望提取的图像表示系数具有优秀的稀疏性与特征保持性，从而只需要融合少量的系数就能获得较优的融合结果。但是，经NSCT变换得到图像低频子带系数的近似为零项十分有限，即不能稀疏的表示图像的低频子带信息，若直接对其融合不利于我们提取源图像的特征。考虑到低频子带包含了图像的主要能量，在很大程度上决定了融合结果的质量，所以我们希望通过提高低频子带系数的稀疏度，以得到更优的融合结果。

发明内容

为了克服现有技术NSCT变化后包含图像主要能量的低频子带系数稀疏度较差，不利于提取有用信息进行融合的不足，本发明提供一种基于NSCT和稀疏表示的图像融合方法，对稀疏度较差的低频子带系数学***自适应调整权重融合；对稀疏度较高的高频方向子带系数采用同一尺度下方向子带绝对值和取最大的方法融合，以捕获源图像中的显著特征，最终提高融合效果。

本发明解决其技术问题所采用的技术方案包括以下步骤：

1训练字典部分：

假设源图像已经经过配准，有K幅大小为M×N的源图像，并分别记作I₁,...,I_K。

(1.1)用NSCT分解每一幅训练图像，经过J级NSCT分解后(J通常为3～5级分解)，得到1个低频子带系数和个高频方向子带系数，其中l_j为尺度j下的方向分解级数。其中，训练图像可以是源图像本身，也可以是与源图像采集方式相同的图像；

(1.2)初始化字典D∈R^n×m，其中，n为字典原子的大小，m为每个子字典的原子数。为了保证字典的过完备性以及计算的复杂度，通常取n＝64,m＝256；

(1.3)对低频子带系数以步长为1，大小为的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵；

(1.4)对上述矩阵用K-SVD算法训练一个字典D，并保存该字典；

2图像融合部分：

(2.1)按照步骤(1.1)的方法用NSCT分解源图像；

(2.2)按照以下6个步骤，融合源图像低频子带系数：

①按照步骤(1.3)中方法将源图像低频子带排列成矩阵V_k,k＝1,...,K；

②将所有源图像的矩阵V_k表示为：

其中，α^C表示共有稀疏表示系数，它包含于所有源图像中；表示第k幅图像的特有稀疏表示系数，它只包含于第k幅源图像中，0表示大小为n×m的全零矩阵，D为步骤(1.4)训练出的字典。

令，

则(1)式可以简化为

V＝D′α  (2)

③为了使得(2)式中α最稀疏，采用正交匹配追踪(OMP)算法求解下式：

$\widehat{\mathrm{\α}}=\arg \underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{\α}\left|\right|}_0,s.t.{\left|\right|D^{\′}\mathrm{\α}-V\left|\right|}_2^2\≤\mathrm{\ϵ}---\left(3\right)$

即，得到α^C和

④按照所有源图像对融合的贡献，将低频子带系数按照下式融合：

${\mathrm{\α}}_f={\mathrm{\α}}^C+\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\ω}}_k{\mathrm{\α}}_k^U,{\mathrm{\ω}}_k=\frac{n_k}{\underset{1\≤i\≤K}{\max }\left(n_k\right)}---\left(4\right)$

其中，n_i表示系数的活动因子，其反应了特征的能量大小，即重要程度：

$n_i={\left|\right|{\mathrm{\α}}_i^U\left|\right|}_1---\left(5\right)$

⑤融合图像低频子带系数可以重构为：

V_f＝Dα_f  (6)

⑥遍历矩阵V_f，将矩阵中的每一列排列成大小的块，再将这些块按照提取的顺序放到的对应位置，并取平均，即对同一位置的子带系数进行累加并除以累加的次数，从而得到源图像低频子带融合系数

(2.3)按照以下2个步骤，融合高频子带系数：

①计算源图像在尺度为2^-l上的方向子带信息，l是在尺度为2^-l分解的方向数，根据经验，2≤l≤4：

$V_l(n,m)=\underset{l\≤i\≤l_j}{\mathrm{\Σ}}\left|V_{l,i}(n,m)\right|---\left(7\right)$

其中，V_l,i(n,m)表示在尺度2^-l、i方向、(n,m)像素位置上的方向子带系数值。

②选取同一尺度下方向子带绝对值和取最大的方法融合：

$V_{l,i}^F(n,m)=V_{l,i}^{k^{*}}(n,m),k^{*}=\underset{l\≤k\≤K}{\arg \max }{V}_l^k(n,m)---\left(8\right)$

其中，1≤l≤J,1≤i≤l_j；和分别表示融合图像或第k^*幅源图像在尺度2^-l、i方向、(n,m)像素位置上的方向子带系数值；k^*表示K幅源图像中，在第l方向子带信息最大的源图像的标号；表示第k幅源图像在第l方向上的子带信息，可根据(7)式定义得到。

(2.4)对融合后的低频子带系数以及融合后的高频子带系数进行NSCT逆变换，得到最终的融合图像F。

本发明的有益效果是：

本发明能有效提高NSCT变换后图像低频子带稀疏度：通过在NSCT域上学习字典，利用稀疏表示求解低频子带的稀疏表示系数，显著提高NSCT变换后图像低频子带稀疏度，更利于提取图像的本质特征，以提高图像融合效果。

本发明较单一的基于NSCT和稀疏表示的融合方法都具有优势：较单一NSCT融合方法的稀疏度更高，更有效提取图像特征，较单一稀疏表示的融合方法具有多尺度、多方向分析能力，更符合人眼观察图像方式。所以，本发明方法较这两大类方法均具有更优的融合效果。

附图说明

图1是提取共有和特有特征图例，其中，(a)是红外源图像，(b)是可见光源图像，(c)是红外图像低频子带，(d)是可见光图像低频子带，(e)是红外图像低频子带特有特征，(f)是可见光图像低频子带特有特征，(g)是红外与可见光低频子带共有特征。

图2是几种方法的融合结果，其中，(a)是DWT方法，(b)是NSCT方法，(c)是SOMP方法，(d)是JSR方法，(e)是本文方法。

图3是本发明图像融合方法流程图。

具体实施方式

一种基于NSCT和稀疏表示的图像融合方法，该方法主要包括训练字典和图像融合两大部分，分述如下：

1训练字典部分：

假设源图像已经经过配准，有K幅大小为M×N的源图像，并分别记作I₁,...,I_K。

(1.1)用NSCT分解每一幅训练图像，经过J级NSCT分解后，得到1个低频子带系数和个高频方向子带系数，其中l_j为尺度j下的方向分解级数。其中，训练图像可以是源图像本身，也可以是与源图像采集方式相同的图像；

(1.2)初始化字典D∈R^n×m，其中，n为字典原子的大小，m为每个子字典的原子数；

(1.3)对低频子带系数以步长为1，大小为的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵；

(1.4)对上述矩阵用K-SVD算法训练一个字典D，并保存该字典；

2图像融合部分：

(2.1)用同步骤(1)中NSCT分解源图像；

(2.2)按照以下5个步骤，融合源图像低频子带系数：

①按照步骤(3)中方法将源图像低频子带排列成矩阵V_k,k＝1,...,K；

②将所有源图像的矩阵V_k表示为：

其中，α^C表示共有稀疏表示系数，它包含于所有源图像中；表示特有稀疏表示系数它只包含于相应一幅源图像中，0表示大小为n×m的全零矩阵。

令，

则(1)式可以简化为

V＝D′α  (2)

③为了使得(2)式中α最稀疏，采用正交匹配追踪(OMP)算法求解下式：

$\widehat{\mathrm{\α}}=\arg \underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{\α}\left|\right|}_0,s.t.{\left|\right|D^{\′}\mathrm{\α}-V\left|\right|}_2^2\≤\mathrm{\ϵ}---\left(3\right)$

即，得到α^C和

④按照所有源图像对融合的贡献，将低频子带系数按照下式融合：

${\mathrm{\α}}_f={\mathrm{\α}}^C+\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\ω}}_k{\mathrm{\α}}_k^U,{\mathrm{\ω}}_k=\frac{n_k}{\underset{1\≤i\≤K}{\max }\left(n_k\right)}---\left(4\right)$

其中，n_i表示系数表示系数的活动因子，其反应了特征的能量大小，即重要程度：

$n_i={\left|\right|{\mathrm{\α}}_i^U\left|\right|}_1---\left(5\right)$

⑤融合图像低频子带系数可以重构为：

V_f＝Dα_f  (6)

⑥遍历矩阵V_f，将矩阵中的每一列排列成大小的块，再将这些块按照提取的顺序放到的对应位置，并取平均，即对同一位置的子带系数进行累加并除以累加的次数，从而得到源图像低频子带融合系数

(2.3)由于NSCT的基具有丰富的方向和形状，能够在高频方向子带上捕获图像中的显著特征，如边缘、线性特征和区域边界。并且这些显著特征在同一尺度所有方向子带上都表现出较大的模值，而其他非显著特征系数模值几乎为零。所以，按照以下2个步骤，融合高频子带系数：

①计算源图像在尺度为2^-l上的方向子带信息：

$V_l(n,m)=\underset{l\≤i\≤l_j}{\mathrm{\Σ}}\left|V_{l,i}(n,m)\right|---\left(7\right)$

其中，V_l,i(n,m)表示在尺度2^-l，i方向，(n,m)像素位置上的方向子带系数值。

②选取同一尺度下方向子带绝对值和取最大的方法融合：

$V_{l,i}^F(n,m)=V_{l,i}^{k^{*}}(n,m),k^{*}=\underset{l\≤k\≤K}{\arg \max }{V}_l^k(n,m)---\left(8\right)$

其中，1≤l≤J,1≤i≤l_j。

(2.4)对融合后的低频子带系数以及融合后的高频子带系数进行NSCT逆变换，得到最终的融合图像F。

下面结合附图和实施例对本发明进一步说明。

实例1.利用稀疏表示提高NSCT低频子带稀疏度并提取源图像特有和共有特征实例

1本实例中字典学习步骤如下：

(1.1)用NSCT分别分解红外与可见光源图像，采用“9-7”塔形分解和“c-d”方向滤波器组，高频层所取的方向数依次为2⁴,2³,2²,2²；

(1.2)初始化字典D∈R^64×256；

(1.3)对低频子带系数以步长为1，大小为8×8的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵，红外低频子带矩阵记为V₁；可见光低频子带矩阵记为V₂；

(1.4)对上述矩阵用K-SVD算法训练一个字典D，K-SVD算法允许误差ε取常用数值0.01，并保存该字典；

2本实例中对低频子带稀疏表示并提取源图像特有和共有特征步骤如下：

(2.1)将红外与可见光源图像的低频子带矩阵表示为：

$\left[\begin{array}{c}{V}_1\\ V_2\end{array}\right]=\left[\begin{array}{ccc}D& D& 0\\ D& 0& D\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}^C\\ {\mathrm{\α}}_1^U\\ {\mathrm{\α}}_2^U\end{array}\right]---\left(9\right)$

其中，α^C表示共有稀疏表示系数，它包含于所有源图像中；和表示特有稀疏表示系数，分别只包含于红外源图像和可见光源图像中，0表示大小为64×256的全零矩阵。

令， $V=\left[\begin{array}{c}{V}_1\\ V_2\end{array}\right],\mathrm{\α}=\left[\begin{array}{c}{\mathrm{\α}}^C\\ {\mathrm{\α}}_1^U\\ {\mathrm{\α}}_2^U\end{array}\right],D^{\′}\left[\begin{array}{ccc}D& D& 0\\ D& 0& D\end{array}\right],$

则(4)式可以简化为

V＝D′α  (10)

(2.2)采用OMP算法求解下式：

$\widehat{\mathrm{\α}}=\arg \underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{\α}\left|\right|}_0,s.t.{\left|\right|D^{\′}\mathrm{\α}-V\left|\right|}_2^2\≤\mathrm{\ϵ}---\left(11\right)$

即，得到α^C、以及

(2.3)分别求得提取的源图像共有特征矩阵VIC以及红外图像特有特征矩阵和可见光图像特有特征矩阵

${\mathrm{VI}}^C=D\×{\mathrm{\α}}^C,{\mathrm{VI}}_1^U=D\×{\mathrm{\α}}_1^U,{\mathrm{VI}}_2^U=D\×{\mathrm{\α}}_2^U---\left(12\right)$

(2.4)分别遍历矩阵和将矩阵中的每一列排列成8×8大小的块，再将这些块按照提取的顺序分别放到I^C，的对应位置，并取平均，即对同一位置的子带系数进行累加并除以累加的次数，从而得到红外与可见光图像共有低频子带系数图像I^C，红外图像和可见光图像各自特有低频子带系数图像和

图1(e)-(g)分别是提取源图像低频子带系数的特有和共有特征后再重构的图像。图像越黑的部分像素值越接近0，零值越多就越稀疏。可以看出，原始红外图像与可见光图像稀疏度较差，没有提取源图像的本质特征。而本发明稀疏度明显提高，并且提取出了红外低频子带中特有的人物和树木轮廓和可见光低频子带中特有的房顶道路等特征，以及他们共有的建筑墙壁和道路等特征，更加利于融合等后续处理。

实例2.本发明图像融合实例

将发明提出的方法与传统基于DWT的图像融合方法和目前性能较为优越的基于NSCT图像融合方法以及基于稀疏表示的图像融合方法SOMP和JSR方法比较。前两种方法是基于变换域的方法，后两种是基于图像域稀疏表示的融合方法。实验中采用240×320大小且经过对准的红外与可见光图像，其中DWT分解的小波类型为3级db4小波，NSCT参数设置与文献相同，即“9-7”塔形分解和“c-d”方向滤波器组，高频层所取的方向数依次为2⁴,2³,2²,2²。稀疏表示的字典大小均为64×256，ε＝0.01。

本实例中基于NSCT和稀疏表示的图像融合方法实现步骤如下：

字典学习步骤与实例1中步骤相同，图像融合部分步骤为：

(1)用NSCT分别分解红外与可见光源图像，NSCT变换参数与学习部分一致，即采用“9-7”塔形分解和“c-d”方向滤波器组，高频层所取的方向数依次为2⁴,2³,2²,2²；

(2)融合源图像低频子带系数：

①通过同实例1中2.1与2.2步骤得到源红外与可见光共有特征系数α^C、源红外图像特有特征系数以及可见光图像特有特征系数

②将低频子带系数按照下式融合：

${\mathrm{\α}}_f={\mathrm{\α}}^C+\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{\mathrm{\ω}}_k{\mathrm{\α}}_k^U,{\mathrm{\ω}}_k=\frac{n_k}{\underset{1\≤i\≤2}{\max }\left(n_k\right)}---\left(13\right)$

其中，n_i表示系数表示系数的活动因子，其反应了特征的能量大小，即重要程度：

$n_i={\left|\right|{\mathrm{\α}}_i^U\left|\right|}_1---\left(14\right)$

⑤融合图像低频子带系数可以重构为：

V_f＝Dα_f  (15)

⑥遍历矩阵V_f，将矩阵中的每一列排列成8×8大小的块，再将这些块按照提取的顺序放到的对应位置，并取平均，即对同一位置的子带系数进行累加并除以累加的次数，从而得到源图像低频子带融合系数

(3)融合高频子带系数：

①求解源图像在尺度为2^-l上的方向子带信息：

$V_l(n,m)=\underset{l\≤i\≤l_j}{\mathrm{\Σ}}\left|V_{l,i}(n,m)\right|---\left(16\right)$

其中，1≤l≤4,1≤i≤l_j(l₁＝2⁴,l₂＝2³,l₃＝2²,l₄＝2²).

②选取同一尺度下方向子带绝对值和取最大的方法融合：

$V_{l,i}^F(n,m)=V_{l,i}^{k^{*}}(n,m),k^{*}=\underset{l\≤k\≤2}{\arg \max }{V}_l^k(n,m)---\left(17\right)$

(4)对融合后的低频子带系数以及融合后的高频子带系数进行NSCT逆变换，得到最终的融合图像F。

融合结果如图2所示，由图2可以看出，DWT融合结果的房屋与周围景物对比度减弱；NSCT融合结果中树木道路层次感较强但是纹理不清晰，人物和房屋也不够突出；SOMP和JSR人物较为突出，但是对树叶等纹理特征明显的区域融合得过于平滑；本文方法则将栏杆、房屋、树叶和道路等地物融合得更为清晰，人物也较为明显，连贯性好，视觉效果最佳。

为了定量评价不同融合方法用于红外和可见光图像融合的性能，本文采用比较方法中指标均方根交叉熵RCE(Rooted Cross Entropy)、Q_W、Q_E和Q_abf进行评价。其中，指标RCE用来评价融合图像与源图像间的综合差异，越小越好；Q_W是源图像与融合图像窗口加权的融合质量评价，Q_E和Q_abf分别从局部和整体反映了融合图像融合源图像边缘的情况，Q₀、Q_W、Q_abf的值均在[0,1]之间，越接近1表明融合质量越好。

表1 几种融合方法的性能指标

表1是几种融合方法的性能指标(其中粗体表示最优的指标值)，观察表1的数据可以看出，相对于在变化域直接融合的方法(DWT和NSCT)与在图像域单尺度基于稀疏表示的方法(SOMP和JSR)，本文提出的方法即能够对源图像进行多尺度分析又能进一步提高图像表示系数的稀疏度，增强了融合图像的细节表现能力，从源图像提取更多的有用信息并加以融合，所以具有更优的融合效果。

Claims (5)

1.一种基于NSCT和稀疏表示的图像融合方法，其特征在于包括下述步骤：

假设源图像已经经过配准，有K幅大小为M×N的源图像，并分别记作I₁,...,I_K；

(1.1)用NSCT分解每一幅训练图像，经过J级NSCT分解后，得到1个低频子带系数和个高频方向子带系数，其中l_j为尺度j下的方向分解级数；

(1.2)初始化字典D∈R^n×m，其中，n为字典原子的大小，m为每个子字典的原子数；

(1.3)对低频子带系数以步长为1，大小为的滑动窗口按照从左上到右下的顺序提取块，再将块拉直并依次排列组成矩阵；

(1.4)对上述矩阵用K-SVD算法训练一个字典D，并保存该字典；

(2.1)按照步骤(1.1)的方法用NSCT分解源图像；

(2.2)按照以下6个步骤，融合源图像低频子带系数：

①按照步骤(1.3)中方法将源图像低频子带排列成矩阵V_k,k＝1,...,K；

②将所有源图像的矩阵V_k表示为：

其中，α^C表示共有稀疏表示系数，它包含于所有源图像中；表示第k幅图像的特有稀疏表示系数，它只包含于第k幅源图像中，0表示大小为n×m的全零矩阵；

令， $V=\left[\begin{array}{c}{V}_1\\ \·\\ \·\\ \·\\ V_K\end{array}\right],$ $\mathrm{\α}=\left[\begin{array}{c}{\mathrm{\α}}^C\\ {\mathrm{\α}}_1^U\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_K^U\end{array}\right],$

则(1)式可以简化为

V＝D′α  (2)

③为了使得(2)式中α最稀疏，采用正交匹配追踪算法求解下式：

$\widehat{\mathrm{\α}}=\arg \underset{\mathrm{\α}}{\min }{\left|\right|\mathrm{\α}\left|\right|}_0,s.t.{\left|\right|D^{\′}\mathrm{\α}-V\left|\right|}_2^2\≤\mathrm{\ϵ}---\left(3\right)$

即，得到α^C和

④按照所有源图像对融合的贡献，将低频子带系数按照下式融合：

${\mathrm{\α}}_f={\mathrm{\α}}^C+\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\ω}}_k{\mathrm{\α}}_k^U,{\mathrm{\ω}}_k=\frac{n_k}{\underset{1\≤i\≤K}{\max }\left(n_k\right)}---\left(4\right)$

其中，n_i表示系数的活动因子，其反应了特征的能量大小，即重要程度：

$n_i={\left|\right|{\mathrm{\α}}_i^U\left|\right|}_1---\left(5\right)$

⑤融合图像低频子带系数可以重构为：

V_f＝Dα_f  (6)

⑥遍历矩阵V_f，将矩阵中的每一列排列成大小的块，再将这些块按照提取的顺序放到的对应位置，并取平均，即对同一位置的子带系数进行累加并除以累加的次数，从而得到源图像低频子带融合系数

(2.3)按照以下2个步骤，融合高频子带系数：

①计算源图像在尺度为2^-l上的方向子带信息，l是在尺度为2^-l分解的方向数：

$V_l(n,m)=\underset{1\≤i\≤l_j}{\mathrm{\Σ}}\left|V_{l,i}(n,m)\right|---\left(7\right)$

其中，V_l,i(n,m)表示在尺度2^-l、i方向、(n,m)像素位置上的方向子带系数值；

②选取同一尺度下方向子带绝对值和取最大的方法融合：

$V_{1,i}^F(n,m)=V_{1,i}^{k^{*}}(n,m),k^{*}=\underset{1\≤k\≤K}{\arg \max }{V}_1^k(n,m)---\left(8\right)$

其中，1≤l≤J,1≤i≤l_j；和分别表示融合图像或第k^*幅源图像在尺度2^-l、i方向、(n,m)像素位置上的方向子带系数值；k^*表示K幅源图像中，在第l方向子带信息最大的源图像的标号；表示第k幅源图像在第l方向上的子带信息；

(2.4)对融合后的低频子带系数以及融合后的高频子带系数进行NSCT逆变换，得到最终的融合图像F。

2.根据权利要求1所述的基于NSCT和稀疏表示的图像融合方法，其特征在于：所述的J取为3～5。

3.根据权利要求1所述的基于NSCT和稀疏表示的图像融合方法，其特征在于：所述的训练图像是源图像本身或是与源图像采集方式相同的图像。

4.根据权利要求1所述的基于NSCT和稀疏表示的图像融合方法，其特征在于：所述的n＝64,m＝256。

5.根据权利要求1所述的基于NSCT和稀疏表示的图像融合方法，其特征在于：所述的l取值范围是2≤l≤4。