CN103049910B - 基于nsct与多通道互信息相结合的图像融合质量评估方法 - Google Patents

Abstract

本发明公开了一种基于NSCT与多通道互信息相结合的图像融合质量评估方法，采用非下采样轮廓波变换方法（NSCT）提取融合前的两幅图像和融合后的图像的层图像，对每层中的多幅图像，计算其对应Contourlet系数的最大值，利用该值对应的图像来计算融合前后每层图像之间的互信息，将每一层图像的权值和互信息相乘，得到融合图像的多模互信息值，借助该值实现融合图像质量评价。本发明使用Contourlet变换与多通道互信息相结合的方法，相较于普通互信息的方法能更精确地反映出融合前与融合后的图像之间的关系，同时与人类视觉感知更加一致，能更客观真实，更真实地出图像融合方法的优劣。

Description

基于NSCT与多通道互信息相结合的图像融合质量评估方法

技术领域

本发明属于融合图像质量评估方法，具体涉及一种基于非下采样轮廓波变换（NSCT）与多通道互信息相结合的医学图像评估方法。

背景技术

图像融合能有效提高图像信息的利用率、提升原始图像的信息量，因此被广泛应用于医学成像领域。图像质量评价对融合的图像质量的评估起到至关重要的作用。图像质量评价能够客观的评价图像融合效果的好坏，能够从客观的角度反映图像融合方法的优劣，从而为判断图像融合算法的优劣提供了客观公正的标准，为图像融合技术的进步起到了积极的推动作用，因此研究有效的图像融合质量评价算法具有重要的理论价值和实际意义。

为了评估不同的融合图像，提出了很多评价算法。比较常用的有互信息（MI）、均方根（RMSE）、通用图像评价指标（UIQI）等。其中，中国专利申请号CN201110432498.5提出了一种基于Gabor特征互信息的全参考图像质量评价方法，该方法利用二维Gabor滤波器提取视觉特征，并通过参考图像和失真图像的互信息来建立图像质量评价函数，取得了符合人类视觉主观认识的评价结果。但是互信息反映的是统计相关性和融合图像中包含源图像的信息量，不能够反映出源图像和融合图像的空间相关性，因而互信息在反应人类视角检测有效性上具有其局限性。而通用图像指数（UIQI）作为一种改进的图像质量评价方法，是通过结合图像的相关损失、亮度失真和对比度失真这三个因素从而建立对于图像失真的模型。尽管这三个因素在图像质量评价中是至关重要的，但该方法不能反映源图像和融合图像的边缘信息。

发明内容

本发明的目的在于提供一种基于NSCT与多通道互信息相结合的图像融合质量评估方法，弥补了仅用互信息无法将图像轮廓的匹配程度纳入评估范围的不足，更真实、客观地反映出融合图像质量的优劣。

基于NSCT与多通道互信息相结合的图像融合质量评估方法，具体为：

（1）采用非下采样轮廓波变换方法提取融合前的两幅图像A,B和融合后的图像F的层图像：

（11）对图像X进行J级非下采样的金字塔滤波得到J个高频率子带图像和一个低频率子带图像，X＝A,B,F；

（12）对图像X的每个高频子带图像进行方向滤波得到张图像，为第个高频率子带图像的方向参数，

（13）将图像X的第个高频率子带图像的张图像上的同一位置像素点(x,y)的contourlet系数取最大值，得到第高频层图像并将低频率子带图像作为低频层图像X₀；

（2）计算融合前的层图像A_j,B_j与融合后的层图像F_j的层间互信息；

（21）在层图像X_j，j＝0,1,…,J中进行窗口逐像素点遍历，窗口每遍历一像素点得到对应的像素点列向量，将层图像A_j与F_j在同一位置像素点的列向量相加得到向量p_i，将层图像B_j与F_j在同一位置像素点的列向量相加得到向量q_i，i＝1,…,N，N为遍历像素点个数；

（22）计算层图像A_j与F_j的协方差矩阵R_AFj以及B_j与F_j的协方差矩阵R_BFj：

$R_{\mathrm{AFj}}=\frac{1}{N}{P}_o{P}_o^T,$ 其中， $P_o=[p_1-\stackrel{\‾}{X},p_2-\stackrel{\‾}{X},p_3-\stackrel{\‾}{X},...,p_N-\stackrel{\‾}{X}],$ $\stackrel{\‾}{X}=\frac{1}{N}{\mathrm{\Σ}}_{i=1}^N{p}_i$

$R_{\mathrm{BFj}}=\frac{1}{N}{Q}_o{Q}_o^T,$ 其中， $Q_o=[q_1-\stackrel{\‾}{Y},q_2-\stackrel{\‾}{Y},q_3-\stackrel{\‾}{Y},...,q_N-\stackrel{\‾}{Y}],$ $\stackrel{\‾}{Y}=\frac{1}{N}{\mathrm{\Σ}}_{i=1}^N{q}_i,$

上标T表示转置；

（23）计算层图像A_j与F_j的互信息RMI(A_j,F_j)以及层图像B_j与F_j的互信息RMI(B_j,F_j)：

$\mathrm{RMI}(A_j,F_j)=\mathrm{Hg}\left(R_{A_j}\right)+\mathrm{Hg}\left(R_{F_{j}1}\right)-\mathrm{Hg}\left(R_{\mathrm{AFj}}\right),$

$\mathrm{RMI}(B_j,F_j)=\mathrm{Hg}\left(R_{B_j}\right)+\mathrm{Hg}\left(R_{F_{j}2}\right)-\mathrm{Hg}\left(R_{\mathrm{BFj}}\right),$

其中，Hg(a)=ln((2πe)^(d/2)det(a)^1/2)， $a=R_{A_j},R_{F_{j}1},R_{F_{j}2},R_{\mathrm{AFj}},R_{B_j},R_{\mathrm{BFj}},$ π为圆周率，e为自然对数的底数，是协方差矩阵R_AFj左上角大小为d/2×d/2的矩阵，R_Fj1是协方差矩阵R_AFj矩阵右下角大小为d/2×d/2的矩阵，是协方差矩阵R_BFj矩阵左上角大小为d/2×d/2的矩阵，R_Fj2是协方差矩阵R_AFj矩阵右下角d/2×d/2的矩阵，d为正整数；

（3）计算融合前的两幅图像A,B和融合后的图像F的多模互信息值MRMI＝MRMI(A F)+MRMI(B，F)，其中 $\mathrm{MRMI}(B,F)=\frac{{\mathrm{\Σ}}_{j=1}^{J+1}{\mathrm{\ω}}_j\×\mathrm{RMI}(B_j,F_j)}{{\mathrm{\Σ}}_{j=1}^{J+1}{\mathrm{\ω}}_j},$ ω_j为第j层图像的权重系数；

（4）以多模互信息值作为融合质量评估值，该评估值越大，说明图像融合质量越高。

所述窗口边长值为奇数。

本发明的技术优点体现在：

利用非下采样contourlet变换对融合前图像和融合后图像进行分解，有利于从多方向和多尺度提取图像的信息，符合人类视觉多通道的特性，而多通道互信息可有效表征多幅图像间的相关性，将非下采样Contourlet变换和多通道互信息相结合，可更真实、更客观地反映出融合图像质量的优劣，为评价各类图像融合算法的性能提供了有效手段。

附图说明

图1是本发明基于非下采样Contourlet变换与多通道互信息相结合的图像评估方法的原理图；

图2是本发明方法详细流程图；

图3是本发明列向量相加示意图；

图4是调制传递函数图；

图5、图6，图7分别是本发明对三组不同的大脑图像的测试结果以及与其它方法的对比示意图，其中，图5a、6a、7a为CT成像图，图5b、6b、7b为MRI成像图，图5c、6c、7c为使用平均方法得到的融合图像，图5d、6d、7d为使用离散小波变换得到的融合图像，5e、6e、7e为使用拉普拉斯金字塔得到的融合图像，5f、6f、7f为使用Ratio金字塔得到的融合图像。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明基于非下采样轮廓波变换（NSCT）与多通道互信息(MRMI)相结合的图像评估方法，包括以下步骤：

（1）提取融合前的两幅图像A,B和融合后的图像F的层图像。

NSCT由两部分构成，分别为非下采样金字塔和方向滤波器，非下采样金字塔的作用为对图像进行滤波，最终依据波段分层；而方向滤波器的作用为对每层图像进行方向滤波得到每个方向的图像。该步骤具体如下：

（11）将尺寸统一为m×n的融合前的图像A(x,y),B(x,y)和融合后的图像F(x,y),经过NSCT的J级非下采样的金字塔滤波，分别得到J个高频率子带图像和一个低频率子带图像，该分层过程为先将图像依据滤波器分解为高频和低频图像，再通过对高频图像进行再一次分解为高频和低频，经历J次分层则可得J个高频率子带图像，一个低频率子带图像。

（12）将每个高频子带图像进行方向滤波，得到（为第个高频率子带图像的方向参数，每层依据需要可设置不同的）张图像。

（13）将第个高频率子带图像的张图像上的同一位置像素点(x,y)的contourlet系数取最大值，设A(x,y),B(x,y)和F(x,y)历经此步骤后得到高频层图像并将低频率子带图像作为低频层图像A₀,B₀,F₀。

（2）计算计算融合前的层图像A_j，B_j与融合后的层图像F_j的层间互信息，j＝0,1,…,J：

（21）选择一个(2r+1)×(2r+1)（正整数0<r<min(m,n),min表示取较小值，m,n为图像尺寸）的窗口，此时窗口的边长为奇数，边长不能为偶数是因为，选取窗口的目的在于获取最中间的像素点与邻域像素点的相互关系，若为偶数，则不存在最中心点。在源图像A_j，B_j以及融合图像F_j上从上到下，从左至右逐像素滑动，当窗口每处于一个像素点位置则可以得到一个尺度为(2r+1)×(2r+1)的列向量，将A_j与F_j在同一位置的列向量相加，得到的向量设为p_i（i＝1,…,N），移动的位置一共为N个点，其中，N＝(m-2r)×(n-2r)，得到矩阵P＝[p₁,p₂,p₃,...,p_N]。将B_j与F_j在同一位置的列向量相加，得到的向量设为q_i（i＝1,…,N），移动的位置一共为N个点，其中，N＝(m-2r)×(n-2r)，得到矩阵Q＝[q₁,q₂,q₃,...,q_N]，如图3所示，方框格子表示在图A和图F中从左上端开始选取（2r+1）(2r+1)的窗口对应点相加（图中r=1），形成尺度为（2r+1）（2r+1）的列向量p₁，然后再向右移动一个像素后同样形成一个尺度为（2r+1）（2r+1）的列向量p₂，以此类推，则一共可移动N次，形成N个列向量，得到矩阵P＝[p₁,p₂,p₃,...,p_N]，对于图B和图F实行此项操作则可得到矩阵Q＝[q₁,q₂,q₃,...,q_N]。

（22）计算矩阵P的协方差矩阵R_AFj:

矩阵Q的协方差矩阵R_BFj:

其中 $P_o=[p_1-\stackrel{\&OverBar;}{X},p_2-\stackrel{\&OverBar;}{X},p_3-\stackrel{\&OverBar;}{X},...,p_N-\stackrel{\&OverBar;}{X}],$ 而其中 $\stackrel{\&OverBar;}{X}=\frac{1}{N}{\mathrm{\&Sigma;}}_{i=1}^N{p}_i$

并且 $Q_o=[q_1-\stackrel{\&OverBar;}{Y},q_2-\stackrel{\&OverBar;}{Y},q_3-\stackrel{\&OverBar;}{Y},...,q_N-\stackrel{\&OverBar;}{Y}],$ 而其中 $\stackrel{\&OverBar;}{Y}=\frac{1}{N}{\mathrm{\&Sigma;}}_{i=1}^N{q}_i$

（23）再根据以下公式得到互信息的值：

$\mathrm{RMI}(A_j,F_j)=\mathrm{Hg}\left(R_{A_j}\right)+\mathrm{Hg}\left(R_{F_{j}1}\right)-\mathrm{Hg}\left(R_{\mathrm{AFj}}\right)$

$\mathrm{RMI}(B_j,F_j)=\mathrm{Hg}\left(R_{B_j}\right)+\mathrm{Hg}\left(R_{F_{j}2}\right)-\mathrm{Hg}\left(R_{\mathrm{BFj}}\right)$

其中Hg(a)=ln((2πe)^(d/2)det(a)^1/2)（a分别等于）

其中R_Aj是R_AFj矩阵左上角大小为d/2×d/2(d为正整数)，R_Fj1是R_AFj矩阵右下角d/2×d/2的矩阵，其中R_Bj是R_BFj矩阵左上角大小为d/2×d/2的矩阵，R_Fj2是R_AFj矩阵右下角d/2×d/2的矩阵（d<min(m,n)），det(a)为计算矩阵a的数值，In为以自然对数e为底的对数。

（3）将根据调制传递函数得到的每一层图片的权值和互信息相乘，得出如下的融合图像的多模互信息值：

MRMI＝MRMI(A,F)+MRMI(B,F)

其中

$\mathrm{MRMI}(A,F)=\frac{{\mathrm{\&Sigma;}}_{j=1}^{J+1}{\mathrm{\&omega;}}_j\&times;\mathrm{RMI}(A_j,F_j)}{{\mathrm{\&Sigma;}}_{j=1}^{J+1}{\mathrm{\&omega;}}_j}$

$\mathrm{MRMI}(B,F)=\frac{{\mathrm{\&Sigma;}}_{j=1}^{J+1}{\mathrm{\&omega;}}_j\&times;\mathrm{RMI}(B_j,F_j)}{{\mathrm{\&Sigma;}}_{j=1}^{J+1}{\mathrm{\&omega;}}_j}$

其中，ω_j表示的第j层图像的权重系数。权重越大表明该层的频率子带对于互信息值的影响就越大。这里的ω_j来自参考文献中的方法，它的基础是调制传递函数（MTF），其图如附录图4所示。当J＝5时，有[ω₁,ω₂,ω₃,ω₄,ω₅,ω₆]=[0.56,0.80,0.92,0.99,0.90,0.40]。

（4）以多模互信息值作为融合质量评估值，该评估值越大，说明图像融合质量越高。

实例：

如图5所示，采用大小为256×256的CT1和MRI1，大小为275×275的CT2和MRI2及CT3和MRI3的图像进行测试，该实例中J=5，该实例中分别使用平均法，离散小波变换法，拉普拉斯金字塔法和比例金字塔法对三组图像分别进行图像融合，都保存为bmp格式，并且用互信息和通用的图像质量指标（UIQI）作为评估标准，对上述融合方法进行评估。

将本发明中所提及的方法以及互信息（MI）和图像质量指标（UIQI）以及作为图像融合评价测度，对用平均法，离散小波变换法，拉普拉斯金字塔法和比例金字塔法融合CT1（图5a）和MRI1（图5b）所得图像如图5中（c，d，e，f）所示进行评估。经测得其各项参数如表1所示：

表1

结合图表可以发现，平均法融合得到的图像信息最丰富，细节多，并且没有过白的现象，而拉普拉斯的方法次之，但也不相上下，而离散小波变换法由于造成了伪影，即使细节保留得好，也不算是一种好的融合方法，而比例金字塔法则是最差的，细节少并且噪声大。三种评价标准的值均为越大，表示融合效果越好，而三种评价方法中只有本方法对于融合方法的评价完全符合上述描述，最符合客观的融合效果评价，并且与视觉效果一致。

将本发明中所提及的方法以及互信息（MI）和图像质量指标（UIQI）以及作为图像融合评价测度，对用平均法，离散小波变换法，拉普拉斯金字塔法和比例金字塔法融合CT2（图6a）和MRI2（图6b）所得图像如图6中（c，d，e，f）所示进行评估。经测得其各项参数如表2所示：

表2

结合图表来看，拉普拉斯金字塔法在该图中的表现最好，不仅轮廓信息完整，内部信息也保留完全，而离散小波变换法则在边缘部分出现波动，稍次之，而平均法保留轮廓信息较完全，但内部细节不够，而比例金子塔法由于噪声的原因，只能处于最末。三种评价标准的值均为越大，表示融合效果越好。本发明方法完全符合上述对于融合方法的好坏排序，拉普拉斯金字塔法获得值最高，离散小波变换法次之，平均法再次，比例金字塔法最末。

将本发明中所提及的方法以及互信息（MI）和图像质量指标（UIQI）以及作为图像融合评价测度，对用对用平均法，离散小波变换法，拉普拉斯金字塔法和比例金字塔法融合CT3（图7a）和MRI3（图7b）所得图像如图7中（c，d，e，f）所示进行评估。经测得其各项参数如表3所示：

表3

如图表所示，平均法获得的信息最丰富，细节完整，表达细腻，而拉普拉斯金字塔法次之，但两者的差距斌不大，而离散小波变换法则又是出现伪影，细节效果不佳，比例金字塔法的噪声掩盖了信息，融合结果。而结合表不难发现本发明提出的方法与视觉效果更为一致，能更真实地评价融合方法的优劣，在对于这两张图的融合结果中，平均法得到的值最高，同时融合效果最好，接下来依次是拉普拉斯金字塔法，离散小波变换法以及比例金字塔法。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.基于NSCT与多通道互信息相结合的图像融合质量评估方法，具体为：

(1)采用非下采样轮廓波变换方法提取融合前的两幅图像A,B和融合后的图像F的层图像：

(11)对图像X进行J级非下采样的金字塔滤波得到J个高频率子带图像和一个低频率子带图像，X＝A,B,F；

(12)对图像X的每个高频子带图像进行方向滤波得到张图像，为第个高频率子带图像的方向参数，

(13)将图像X的第个高频率子带图像的张图像上的同一位置像素点(x,y)的contourlet系数取最大值，得到第高频层图像并将低频率子带图像作为低频层图像X₀；

(2)计算融合前的层图像A_j,B_j与融合后的层图像F_j的层间互信息；

(21)在层图像X_j，j＝0,1,…,J中进行窗口逐像素点遍历，窗口每遍历一像素点得到对应的像素点列向量，将层图像A_j与F_j在同一位置像素点的列向量相加得到向量p_i，将层图像B_j与F_j在同一位置像素点的列向量相加得到向量q_i，i＝1,…,N，N为遍历像素点个数；

(22)计算层图像A_j与F_j的协方差矩阵R_AFj以及B_j与F_j的协方差矩阵R_BFj：

$R_{\mathrm{AFj}}=\frac{1}{N}{P}_o{P}_o^T,$ 其中， $P_o=[p_1-\stackrel{\&OverBar;}{X},p_2-\stackrel{\&OverBar;}{X},p_3-\stackrel{\&OverBar;}{X},...,p_N-\stackrel{\&OverBar;}{X}],\stackrel{\&OverBar;}{X}=\frac{1}{N}{\mathrm{\&Sigma;}}_{i=1}^N{p}_i$

$R_{\mathrm{BFj}}=\frac{1}{N}{Q}_o{Q}_o^T,$ 其中， $Q_o=[q_1-\stackrel{\&OverBar;}{Y},q_2-\stackrel{\&OverBar;}{Y},q_3-\stackrel{\&OverBar;}{Y},...,q_N-\stackrel{\&OverBar;}{Y}],\stackrel{\&OverBar;}{Y}=\frac{1}{N}{\mathrm{\&Sigma;}}_{i=1}^N{q}_i,$

上标T表示转置；

(23)计算层图像A_j与F_j的互信息RMI(A_j,F_j)以及层图像B_j与F_j的互信息RMI(B_j,F_j)：

$\mathrm{RMI}(A_j,F_j)=\mathrm{Hg}\left(R_{A_j}\right)+\mathrm{Hg}\left(R_{F_{j}1}\right)-\mathrm{Hg}\left(R_{\mathrm{AFj}}\right),$

$\mathrm{RMI}(B_j,F_j)=\mathrm{Hg}\left(R_{B_j}\right)+\mathrm{Hg}\left(R_{F_{j}2}\right)-\mathrm{Hg}\left(R_{\mathrm{BFj}}\right),$

其中，Hg(a)＝ln((2πe)^(d/2)det(a)^1/2)， $a=R_{A_j},R_{F_{j}1},R_{F_{j}2},R_{\mathrm{AFj}},R_{B_j},R_{\mathrm{BFj}},$ π为圆周率，e为自然对数的底数，是协方差矩阵左上角大小为d/2×d/2的矩阵，是协方差矩阵R_AFj矩阵右下角大小为d/2×d/2的矩阵，是协方差矩阵R_BFj矩阵左上角大小为d/2×d/2的矩阵，是协方差矩阵R_AFj矩阵右下角d/2×d/2的矩阵，d为正整数，det(a)为计算矩阵a的行列式；

(3)计算融合前的两幅图像A,B和融合后的图像F的多模互信息值MRMI＝MRMI(A,F)+MRMI(B,F)，其中 $\mathrm{MRMI}(B,F)=\frac{{\mathrm{\&Sigma;}}_{j=0}^J{\mathrm{\&omega;}}_j\&times;\mathrm{RMI}(B_j,F_j)}{{\mathrm{\&Sigma;}}_{j=0}^J{\mathrm{\&omega;}}_j},$ ω_j为第j层图像的权重系数；

(4)以多模互信息值作为融合质量评估值，该评估值越大，说明图像融合质量越高。

2.根据权利要求1所述的图像融合质量评估方法，其特征在于，所述窗口边长值为奇数。