图像融合质量评价方法MSSIM、MS-SSIM、FS、Qmi、Qabf与VIFF（三）

1 前言

在前两篇文章中介绍的是比较经典比较常见的融合图像质量评价指标，可以点击图像融合质量评价方法SSIM、PSNR、EN、MSE与NRMSE（一）和图像融合质量评价方法AG、SF、STD、MI与NMI（二）查看。后续更新了图像融合质量评价方法FMI（四）
下面介绍一下，比较新颖的评价指标。

2 融合评价指标介绍

2.1 结构相似性（SSIM）相关

首先就是MSSIM，这个很容易理解就是源图像A与融合图像F的SSIM与源图像B与融合图像F的SSIM的平均值，也有写成$SSI{M}_a(F)$,公式如下：$$MSSIM=(SSIM(A,F)+SSIM(B,F))$$其次是多层级结构相似性(Multiscale structural similarity),简写为MS-SSIM。MS-SSIM能更好地与人眼视觉系统的视觉感知相一致，并且在一定的尺度下，评价效果优于SSIM，公式如下：$$MSSSIM(Z,K)={\left[l_M(Z,K)\right]}^{{\alpha }_M}\prod _{i=1}^M{\left[s_i(Z,K)\right]}^{{\beta }_i}{\left[z_i(Z,K)\right]}^{{\gamma }_i}$$具体可见Multiscale structural similarity for image quality assessment

2.2 互信息（MI）相关

根据文章Information measure for performance of image fusion中指出，用于评估融合图像的互信息为：$$MI=I_{AF}+I_{BF}$$$I_{AF}$与$I_{BF}$分别表示源图像A与融合图像F以及源图像B与融合图像F的互信息。
再根据IEEE公布2002年的论文Fusion performance measures and a lifting wavelet transform based algorithm for image fusion引入基于互信息（MI）的融合因子（Fusion Factor）与融合对称因子（Fusion Symmetry）的概念。其中融合因子FF为：$$FF=I_{AF}+I_{BF}$$较高的FF值表示融合图像包含在两个图像中都存在的相当大量的信息。但是，高FF值并不意味着来自两个图像的信息都是对称融合的。所以引入融合对称因子FS：$$FS=abs(\frac{I_{AF}}{I_{AF}+I_{BF}}-0.5)$$FS的值越小表示融合算法的效果越好。
基于Image fusion based on multi-scale guided filters这篇文章中提出了归一化互信息$Q_{MI}$,公式如下：$$Q_{MI}=2\left[\frac{I_{AF}}{H_A+H_F}+\frac{I_{BF}}{H_B+H_F}\right]$$$Q_{MI}$的值越大表示从源图像中获取的信息越多，融合的效果越好。

2.3 基于视觉信息保真度的指标VIFF

VIFF（The visual information fidelity for fusion）是基于视觉信息保真度提出的衡量融合图像质量的指标，由于太复杂了，详见A new image fusion performance metric based on visual information fidelity。VIFF的值越大，表示融合图像质量越好。

2.4 融合质量Qabf

Qabf是一种新颖的融合图像客观非参考质量评估指标，得到Qabf的算法利用局部度量来估计来自输入的显着信息在融合图像中的表现程度，Qabf的值越高，表示融合图像的质量越好。计算公式如下$$Q(a,b,f)=\frac{1}{\left|W\right|}\sum _{\omega \in W}(\lambda (\omega )Q_0(a,f|\omega )+(1-\lambda (\omega ))Q_0(b,f|\omega ))$$详细可见见论文A new quality metric for image fusion

3 代码实现

3.1 多层次结构相似性MS-SSIM matlab实现

来源网上代码，并不会写matlab，但是勉强各种参考可以用：

function overall_mssim = msssim(img1, img2, K, win, level, weight, method)

% Multi-scale Structural Similarity Index (MS-SSIM)
% Z. Wang, E. P. Simoncelli and A. C. Bovik, "Multi-scale structural similarity
% for image quality assessment," Invited Paper, IEEE Asilomar Conference on
% Signals, Systems and Computers, Nov. 2003

if (nargin < 2 || nargin > 7)
   overall_mssim = -Inf;
   return;
end

if (~exist('K'))
   K = [0.01 0.03];
end

if (~exist('win'))
   win = fspecial('gaussian', 11, 1.5);
end

if (~exist('level'))
   level = 5;
end

if (~exist('weight'))
   weight = [0.0448 0.2856 0.3001 0.2363 0.1333];
end

if (~exist('method'))
   method = 'product';
end

if (size(img1) ~= size(img2))
   overall_mssim = -Inf;
   return;
end

[M N] = size(img1);
if ((M < 11) || (N < 11))
   overall_mssim = -Inf;
   return
end

if (length(K) ~= 2)
   overall_mssim = -Inf;
   return;
end

if (K(1) < 0 || K(2) < 0)
   overall_mssim = -Inf;
   return;
end
  
[H W] = size(win);

if ((H*W)<4 || (H>M) || (W>N))
   overall_mssim = -Inf;
   return;
end
   
if (level < 1)
   overall_mssim = -Inf;
   return
end


min_img_width = min(M, N)/(2^(level-1));
max_win_width = max(H, W);
if (min_img_width < max_win_width)
   overall_mssim = -Inf;
   return;
end

if (length(weight) ~= level || sum(weight) == 0)
   overall_mssim = -Inf;
   return;
end

if (method ~= 'wtd_sum' & method ~= 'product')
   overall_mssim = -Inf;
   return;
end

downsample_filter = ones(2)./4;
im1 = double(img1);
im2 = double(img2);
for l = 1:level
   [mssim_array(l) ssim_map_array{l} mcs_array(l) cs_map_array{l}] = ssim_index_new(im1, im2, K, win);
%   [M N] = size(im1);
   filtered_im1 = imfilter(im1, downsample_filter, 'symmetric', 'same');
   filtered_im2 = imfilter(im2, downsample_filter, 'symmetric', 'same');
   clear im1, im2;
   im1 = filtered_im1(1:2:end, 1:2:end);
   im2 = filtered_im2(1:2:end, 1:2:end);
end

if (method == 'product')
%   overall_mssim = prod(mssim_array.^weight);
   overall_mssim = prod(mcs_array(1:level-1).^weight(1:level-1))*(mssim_array(level).^weight(level));
else
   weight = weight./sum(weight);
   overall_mssim = sum(mcs_array(1:level-1).*weight(1:level-1)) + mssim_array(level).*weight(level);
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

这个是调用的文件：

function Q = assess(im1,im2,fused)
    im1=imread('C:\Users\DELL\Desktop\code\1.png')
    im2=imread('C:\Users\DELL\Desktop\code\2.png')
    fused=imread('C:\Users\DELL\Desktop\code\fu.png')
    im1=double(im1);
    im2=double(im2);
    fused = double(fused);
    Q = [];

    Q(1)=msssim(im1,fused);
end
1
2
3
4
5
6
7
8
9
10
11

3.2 互信息MI相关指标 python实现

根据上面所给的公式，很容易可以得到代码：

import cv2
import numpy as np
import math

from MI import mutualInfo
from skimage.measure import shannon_entropy

# 融合因子
def fusionFactor(img1,img2,fu_img):
    
    ff = mutualInfo(img1,fu_img)+mutualInfo(img2,fu_img)

    return ff

# 融合对称因子
def fusionSymmetry(img1,img2,fu_img):

	MI1 = mutualInfo(img1,fu_img)
	MI2 = mutualInfo(img2,fu_img)

	fs = math.abs(MI1/(MI1+MI2) - 0.5)

	return fs

def QMI(img1,img2,fu_img):

    en1 = shannon_entropy(img1)
    en2 = shannon_entropy(img2)
    en_fu = shannon_entropy(fu_img)

    mi1 = mutualInfo(img1,fu_img)
    mi2 = mutualInfo(img2,fu_img)
    
    Qmi = 2.0*((mi1/(en1+en_fu))+(mi2/(en2+en_fu)))

    return Qmi


if __name__ == '__main__':
    img1 = cv2.imread('1.png',0)
    img2 = cv2.imread('2.png',0)
    img_fu = cv2.imread('3.png',0)

    print(fusionFactor(img1,img2,img_fu))
    print(fusionSymmetry(img1,img2,img_fu))
    print(QMI(img1,img2,img_fu))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

3.5 VIFF matlab实现

在github上面找到HarrisXia/image-fusion-evaluation中的VIFF_Public.m文件，如下：

%% -----------COPYRIGHT NOTICE STARTS WITH THIS LINE------------
% Copyright (c) 2012, Yu Han,(Chinese Name: HAN Yu) All rights reserved.
%       The name of this code is "image fusion performance metric based on visual information fidelity".
% Permission to use and copy this software and its documentation for educational and 
% research purposes only and without fee is hereby granted, provided that this 
% copyright notice and the original authors names appear on all copies and supporting 
% documentation. 
%   The authors are acknowledged in any publication that reports research using this software.
%   The work is to be cited in the bibliography as:
%   	[]Yu Han, Yunze Cai, Yin Cao, Xiaoming Xu, A new image fusion performance metric 
%   	based on visual information fidelity, information fusion, Volume 14, Issue 2, April 2013, Pages 127–135
%   This code shall not be used, rewritten, or adapted as the basis of a commercial 
% software or hardware product without hand-writing permission of the authors. The authors 
% make no representations about the suitability of this software for any purpose. It is 
% provided "as is" without express or implied warranty.
%% -----------COPYRIGHT NOTICE ENDS WITH THIS LINE------------
function F=VIFF_Public(Im1,Im2,ImF)
% Cite this work as
% []Yu Han, Yunze Cai, Yin Cao, Xiaoming Xu, A new image fusion performance metric based on visual information fidelity, information fusion, Volume 14, Issue 2, April 2013, Pages 127–135
% input:
%       Im1, source image 1
%       Im2, source image 2
%       ImF, fused image
% output:
%       F, fusion assessment value
%
% visual noise
sq=0.005*255*255;
% error comaprison parameter
C=1e-7;

[r,s,l]=size(Im1);
%color space transformation
if l==3
    cform = makecform('srgb2lab');
    T1 = applycform(Im1,cform);
    T2 = applycform(Im2,cform);
    TF = applycform(ImF,cform);
    Ix1=T1(:,:,1);
    Ix2=T2(:,:,1);
    IxF=TF(:,:,1); 
else
    Ix1=Im1;
    Ix2=Im2;
    IxF=ImF;
end

T1p=double(Ix1);
T2p=double(Ix2);
Trp=double(IxF);

p=[1,0,0.15,1]./2.15;
[T1N,T1D,T1G]=ComVidVindG(T1p,Trp,sq);
[T2N,T2D,T2G]=ComVidVindG(T2p,Trp,sq);
VID=[];
VIND=[];
%i multiscale image level
for i=1:4
    M_Z1=cell2mat(T1N(i));
    M_Z2=cell2mat(T2N(i));
    M_M1=cell2mat(T1D(i));
    M_M2=cell2mat(T2D(i));
    M_G1=cell2mat(T1G(i));
    M_G2=cell2mat(T2G(i));
    L=M_G1<M_G2;
    M_G=M_G2;
    M_G(L)=M_G1(L);
    M_Z12=M_Z2;
    M_Z12(L)=M_Z1(L);
    M_M12=M_M2;
    M_M12(L)=M_M1(L);
    
    VID=sum(sum((M_Z12+C)));
    VIND=sum(sum((M_M12+C)));
    F(i)=VID/VIND;
end
F=sum(F.*p);


function [Tg1,Tg2,Tg3]=ComVidVindG(ref,dist,sq)
% this part is mainly from the work:
% [] H.R.Sheikh and A.C.Bovik, Image information and visual quality[J], IEEE Transactions on Image Processing 15(2), pp. 430–444, 2006.
% And we have a little revision in our code
% input:
%       ref, source image
%       dist,fused image
%       sq, visual noise
% output:
%       Tg1, the matrix of visual information with distortion information (VID)
%       Tg2, the matrix of visual information without distortion information (VIND)
%       Tg3, the matrix of scalar value gi
sigma_nsq=sq;

for scale=1:4
   
    N=2^(4-scale+1)+1;
    win=fspecial('gaussian',N,N/5);
    
    if (scale >1)
        ref=filter2(win,ref,'valid');
        dist=filter2(win,dist,'valid');
        ref=ref(1:2:end,1:2:end);
        dist=dist(1:2:end,1:2:end);
    end
    
    mu1   = filter2(win, ref, 'valid');
    mu2   = filter2(win, dist, 'valid');
    mu1_sq = mu1.*mu1;
    mu2_sq = mu2.*mu2;
    mu1_mu2 = mu1.*mu2;
    sigma1_sq = filter2(win, ref.*ref, 'valid') - mu1_sq;
    sigma2_sq = filter2(win, dist.*dist, 'valid') - mu2_sq;
    sigma12 = filter2(win, ref.*dist, 'valid') - mu1_mu2;
    
    sigma1_sq(sigma1_sq<0)=0;
    sigma2_sq(sigma2_sq<0)=0;
    
    g=sigma12./(sigma1_sq+1e-10);
    sv_sq=sigma2_sq-g.*sigma12;
    
    g(sigma1_sq<1e-10)=0;
    sv_sq(sigma1_sq<1e-10)=sigma2_sq(sigma1_sq<1e-10);
    sigma1_sq(sigma1_sq<1e-10)=0;
    
    g(sigma2_sq<1e-10)=0;
    sv_sq(sigma2_sq<1e-10)=0;
    
    sv_sq(g<0)=sigma2_sq(g<0);
    g(g<0)=0;
    sv_sq(sv_sq<=1e-10)=1e-10;
    
     G(scale)={g};
     VID=log10(1+g.^2.*sigma1_sq./(sv_sq+sigma_nsq));
     VIND=log10(1+sigma1_sq./sigma_nsq);
     Num(scale)={VID};
     Den(scale)={VIND};    
end
Tg1=Num;
Tg2=Den;
Tg3=G;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

在matlab中写一个调用函数的文件，即可使用，最好不要单独下载使用。

3.4 Qabf matlab实现

同样来源网上，matlab实现，将图片地址改掉后直接运行即可：

%    Objective Image Fusion Performance Measure 
%    C.S.Xydeas and V.Petrovic 
%    Electronics Letters 
 
%    edit by Richang Hong; May 27 07 
 
function output = Qabf(strA, strB, strF) 
% strA and strB is the source images and strF is the fusion result 
 
% model parameters 
    L=1; Tg=0.9994;kg=-15;Dg=0.5;Ta=0.9879;ka=-22;Da=0.8;     
     
% Sobel Operator 
h1=[1 2 1;0 0 0;-1 -2 -1]; h2=[0 1 2;-1 0 1;-2 -1 0]; h3=[-1 0 1;-2 0 2;-1 0 1]; 
% if y is the response to h1 and x is the response to h3; 
% then the intensity is sqrt(x^2+y^2) and orientation is arctan(y/x); 
strA=imread('C:\Users\DELL\Desktop\code\1.png')
strB=imread('C:\Users\DELL\Desktop\code\2.png')
strF=imread('C:\Users\DELL\Desktop\code\fu.png')
%strF=imread('D:\hb-123456\深度学习U-NET的研究\训练\许-203融合结果\case10_fusion.tif')
pA = double(strA);
pB = double(strB); 
pF = double(strF);
%pA = imread(strA); pA = double(pA); 
%pB = imread(strB); pB = double(pB); 
%pF = imread(strF); pF = double(pF); 
 
SAx = conv2(pA,h3,'same'); SAy = conv2(pA,h1,'same'); 
gA = sqrt(SAx.^2 + SAy.^2);  
[M,N] = size(SAx); aA = zeros(M,N); 
for i=1:M 
    for j=1:N 
        if ( SAx(i,j) == 0 ) aA(i,j) = pi/2; 
        else 
            aA(i,j) = atan(SAy(i,j)/SAx(i,j)); 
        end 
    end 
end 
 
SBx = conv2(pB,h3,'same'); SBy = conv2(pB,h1,'same'); 
gB = sqrt(SBx.^2 + SBy.^2);  
[M,N] = size(SBx); aB = zeros(M,N); 
for i=1:M 
    for j=1:N 
        if ( SBx(i,j) == 0 ) aB(i,j) = pi/2; 
        else 
            aB(i,j) = atan(SBy(i,j)/SBx(i,j)); 
        end 
    end 
end 
 
SFx = conv2(pF,h3,'same'); SFy = conv2(pF,h1,'same'); 
gF = sqrt(SFx.^2 + SFy.^2);  
[M,N] = size(SAx); aF = zeros(M,N); 
for i=1:M 
    for j=1:N 
        if ( SFx(i,j) == 0 ) aF(i,j) = pi/2; 
        else 
            aF(i,j) = atan(SFy(i,j)/SFx(i,j)); 
        end 
    end 
end 
 
% the relative strength and orientation value of GAF,GBF and AAF,ABF; 
GAF = zeros(M,N); AAF = zeros(M,N); QgAF = zeros(M,N); QaAF = zeros(M,N); QAF = zeros(M,N); 
for i=1:M 
    for j=1:N 
         
        if ( gA(i,j) > gF(i,j))  GAF(i,j) = gF(i,j)/gA(i,j); 
        else 
            if ( gA(i,j) == gF(i,j) )  GAF(i,j) = gF(i,j); 
            else 
                GAF(i,j) = gA(i,j) / gF(i,j); 
            end 
        end  
        AAF(i,j) = 1 - abs(aA(i,j)-aF(i,j))/(pi/2); 
         
        QgAF(i,j) = Tg / (1 + exp(kg*( GAF(i,j) - Dg ))); 
        QaAF(i,j) = Ta / (1 + exp(ka*( AAF(i,j) - Da ))); 
         
        QAF(i,j) = QgAF(i,j) * QaAF(i,j); 
    end 
end 
 
GBF = zeros(M,N); ABF = zeros(M,N); QgBF = zeros(M,N); QaBF = zeros(M,N); QBF = zeros(M,N); 
for i=1:M 
    for j=1:N 
         
        if ( gB(i,j) > gF(i,j))  GBF(i,j) = gF(i,j)/gB(i,j); 
        else 
            if ( gB(i,j) == gF(i,j) )  GBF(i,j) = gF(i,j); 
            else 
                GBF(i,j) = gB(i,j) / gF(i,j); 
            end 
        end  
        ABF(i,j) = 1 - abs(aB(i,j)-aF(i,j))/(pi/2); 
         
        QgBF(i,j) = Tg / (1 + exp(kg*( GBF(i,j) - Dg ))); 
        QaBF(i,j) = Ta / (1 + exp(ka*( ABF(i,j) - Da ))); 
         
        QBF(i,j) = QgBF(i,j) * QaBF(i,j); 
    end 
end 
 
% compute the QABF 
deno = sum(sum( gA + gB )); 
nume = sum(sum( QAF.*gA + QBF.*gB )); 
output = nume / deno; 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

4 总结

目前更新到这里，如果后面还收集了更多指标和代码实现，会及时更新的…