数字图像融合相关评价指标

数字图像融合相关评价指标

数字图像融合质量评价指标：标准差、平均梯度、边缘强度、空间频率、信息熵、交叉熵、互信息、结构相似度等。
-----------------这里只罗列部分代码哦-----------------------
有关于详细论文如下：融合图像质量评价指标的相关性分析及性能评估–张小利

1、标准差：标准差用于表述像素灰度值的离散程度，可以反应图像细节信息的丰富程度，值越大表示，细节越丰富。
2、平均梯度：平均梯度用于表征图像的清晰程度，计算的结果值越大，表明融合图像越清晰。

//这里是MATLAB代码！
/*
author:Anter
*/
function outval = avg_gradient(img) 
% OUTVAL = AVG_GRADIENT(IMG) 
%  img=im2uint8(img);
if nargin == 1 
    img = double(img); 
    % Get the size of img 
    [r,c,b] = size(img); 
     
    dx = 1; 
    dy = 1; 
    for k = 1 : b 
        band = img(:,:,k); 
        [dzdx,dzdy] = gradient(band,dx,dy); 
        s = sqrt((dzdx .^ 2 + dzdy .^2) ./ 2); 
        g(k) = sum(sum(s)) / ((r - 1) * (c - 1)); 
    end 
    outval = mean(g); 
else 
    error('Wrong number of input!'); 
end
end

3、边缘强度：

//这里是MATLAB代码！
/*
author:Anter
*/
function outval = edge_intensity(img) 
if nargin == 1 
    img = double(img); 
    % Create horizontal sobel matrix 
    w = fspecial('sobel'); 
     
    % Get the size of img 
    [r c k] = size(img); 
     
    gx = imfilter(img,w,'replicate'); 
    gy = imfilter(img,w','replicate'); 
     
    for m = 1 : r 
        for n = 1 : c 
            for q = 1 : k 
                g(m,n,q) = sqrt(gx(m,n,q)*gx(m,n,q) + gy(m,n,q)*gy(m,n,q)); 
            end 
        end 
    end 
    outval = mean(mean(mean(g))); 
 else 
    error('Wrong number of input!'); 
 end
end

4、峰值信噪比：
峰值信噪比是用来衡量结果图中有效信息和噪声比例的变量，可以定量的表述融合图像的失真程度。

//这里是MATLAB代码！
/*
author:Anter

*/
function psnr=PSNR(A,F)

if size(A,3)==3
    A=rgb2gray(A);
end
if size(F,3)==3
    F=rgb2gray(F);
end
[height width]=size(A);
A = im2double(A);
F = im2double(F);

sigma1=0;
for i=1:height
    for j=1:width
        sigma1=sigma1+(A(i,j)-F(i,j))^2;
    end
end

MSE=(sigma1/(height*width));  %均方误差
psnr=10*log10((255^2)/MSE)

4、空间频率
空间频率主要反应图像灰度变化率，是图像边缘信息量的评判指标。

//这里是MATLAB代码！
function SF=space_frequency(X)
% X=im2uint8(X);
X=double(X);
[n0,n1]=size(X);%%%%   X是一个方阵
X=double(X);                          %空间频率;
RF=0;
CF=0;

for fi=1:n0
    for fj=2:n1
        RF=RF+(X(fi,fj)-X(fi,fj-1)).^2;
    end
end

RF=RF/(n0*n1);

for fj=1:n1
    for fi=2:n0
        CF=CF+(X(fi,fj)-X(fi-1,fj)).^2;
    end
end

CF=CF/(n0*n1);%%%%可以思考，空间频率是不是只描述一个方阵

SF=sqrt(RF+CF)
end

5、信息熵
信息熵用于衡量图像中信息的丰富程度，属于正向指标，值越大，表示融合图像包含的信息越丰富。

//这里是MATLAB代码！

function S = shannon(img) 

% img=im2uint8(img);
I=img;
%I=double(I);
[C,R]=size(I);      %求图像的规格
Img_size=C*R;       %图像像素点的总个数
L=256;              %图像的灰度级
H_img=0;
nk=zeros(L,1);
for i=1:C
    for j=1:R
        Img_level=I(i,j)+1;                 %获取图像的灰度级           
        nk(Img_level)=nk(Img_level)+1;      %统计每个灰度级像素的点数
    end
end
for k=1:L
    Ps(k)=nk(k)/Img_size;                  %计算每一个灰度级像素点所占的概率
    if Ps(k)~=0;                           %去掉概率为0的像素点
    H_img=-Ps(k)*log2(Ps(k))+H_img;        %求熵值的公式
    S=H_img;
    end
end
end

6、结构相似度(SSIM)
结构相似性用于表述结果图与源图像结构的相似程度，主要由亮度、对比度损失和相关性损失三部分构成，融合图像对源图像的结构保存越多，结构相似性越大。

//这里是MATLAB代码！
if size(img1,3)==3
    img1=rgb2gray(img1);
end
if size(img2,3)==3
    img2=rgb2gray(img2);
end
[M N] = size(img1);

img1 = double(img1);
img2 = double(img2);
K = [0.01 0.03]; window = fspecial('gaussian', 11, 1.5);L = 255;
% automatic downsampling
f = max(1,round(min(M,N)/256));
%downsampling by f
%use a simple low-pass filter 
if(f>1)
    lpf = ones(f,f);
    lpf = lpf/sum(lpf(:));
    img1 = imfilter(img1,lpf,'symmetric','same');
    img2 = imfilter(img2,lpf,'symmetric','same');

    img1 = img1(1:f:end,1:f:end);
    img2 = img2(1:f:end,1:f:end);
end

C1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
window = window/sum(sum(window));

mu1   = filter2(window, img1, 'valid');
mu2   = filter2(window, img2, 'valid');
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;

if (C1 > 0 & C2 > 0)
   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
else
   numerator1 = 2*mu1_mu2 + C1;
   numerator2 = 2*sigma12 + C2;
	denominator1 = mu1_sq + mu2_sq + C1;
   denominator2 = sigma1_sq + sigma2_sq + C2;
   ssim_map = ones(size(mu1));
   index = (denominator1.*denominator2 > 0);
   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
   index = (denominator1 ~= 0) & (denominator2 == 0);
   ssim_map(index) = numerator1(index)./denominator1(index);
end
SSIM = mean2(ssim_map)

7、边缘信息评价因子
Q用于衡量融合图像中有多少边缘信息量来源于源图像，其值越接近1，表示来源于源图像的边缘信息转移就越多。

function Qabf =QABF(pA,pB,pF)
% 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);
if size(pA,3)==3
    pA=rgb2gray(pA);
end
if size(pB,3)==3
    pB=rgb2gray(pB);
end
if size(pF,3)==3
    pF=rgb2gray(pF);
end
pA = im2double(pA);
pB = im2double(pB);
pF = im2double(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 ));
Qabf = nume / deno

----------以上是图像融合部分评价指标，亲测有效哦！----------------