Objective Evaluation Index of image
图像质量客观评价指标
在做红外图像细节增强算法研究时，很重要一点就是要对经过算法处理的图像结果进行评价，分成两种评价方法。一种是视觉效果评价：主观的人眼观察，主要是通过观察者能否看到更多图像细节，给人的感觉是否真实自然，与其他算法处理的图像的亮度和对比度的比较。在论文中的呈现方式，就是把各个算法的处理结果和原图一起呈现给观察者，再对图像某一具体区域进行详细分析，如图像中的文字，黑暗区域，栏杆树木和建筑物等细节信息进行比较，看那种算法的结果更清晰，暴露出更多的细节。当然咯，主观评价方式往往取决于观察者，并不是一种客观可靠的方式，往往需要几种公认的图像质量客观评价指标，对红外图像细节增强算法的处理结果进行比较。通过阅读大量的SCI文献，如IEEE Transactions on Image Processing，Optic Engineering, Applied Optics, Optik, Infrared Physics & Technology等多种权威期刊的图像处理论文精读，本文选取了四种经典的客观评价指标：RMSC, SSIM, FSIM, MOS
建议阅读另一篇博客，关于matlab批量输出，会更透彻，链接如下：
https://www.cnblogs.com/Xwangzi66/p/15171728.html
一 RMSC
The Root-Mean-Square Contrast (RMSC) index，即均方根对比度，用来衡量图像的对比度，其值越大对比度越高。其中， 为图像总像素的均值。I为图像某点的像素。M和N分别为图像的高和宽。
The MATLAB code of RMSC is following:
直接根据定义编写matlab代码即可，但是为了更实用，避免手动一张张图像的仿真，很耗时！我们可以把它封装成一个函数，并且推广多张图像一次显示RMSC的值。思路：对于批量输出的算法，一般是把算法封装成一个函数，同时结合循环，直接调用即可。
RMSC的计算代码分成两部分：一个是计算RMSC值的功能函数。另一个是顶层，调用功能函数，批量输出多张图像的RMSC的值。
一个是功能函数，计算RMSC：
function rmsc_out = RMSC_cal(I) %load('S1.mat'); %I = double(S); % I = imread('./DATA/Ablation_wFilter.tif'); % I = imread('./DATA/LYF2.jpg'); %I= rgb2gray(I);%将 RGB 图像或颜色图转换为灰度图 I= double(I); avg=mean2(I); %求图像均值 [m,n]=size(I); s=0; for x=1:m for y=1:n s=s+(I(x,y)-avg)^2; %求得所有像素与均值的平方和。 end end %故RMSC的MATLAB代码如下： rmsc_out=sqrt(s/(m*n)); %利用方差公式求得. sqrt返回数组 X 的每个元素的平方根 end
上面为计算RMSC的值，下面为顶层函数，直接输出五张图像的RMSC值，且得到五张图像的平均值，同理可以推广到成百上千张图像，代码如下
1 clc; 2 close all; 3 clear; 4 %--------without the ROGABF 5 I1=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF1.png'); 6 I2=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF2.png'); 7 I3=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF3.png'); 8 I4=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF4.png'); 9 %324*256 10 %512*512 11 I5=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF11.png'); 12 %--------without the ROGABF 13 %--------without the wFilter 14 % I1=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter1.png'); 15 % I2=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter2.png'); 16 % I3=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter3.png'); 17 % I4=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter4.png'); 18 % %324*256 19 % %512*512 20 % I5=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter11.png'); 21 %------without the wFilter 22 %output five result of RMSC at a time 23 RMSC_out = zeros(1,5); 24 I = cell(1,5); 25 I{1,1}=I1;I{1,2}=I2; 26 I{1,3}=I3;I{1,4}=I4;I{1,5}=I5; 27 for k = 1:5 28 temp = I{k}; 29 RMSC_out(k) = RMSC_cal(temp); 30 end 31 RMSC_out; 32 mean(RMSC_out) 33 %output five result of RMSC at a time
二 SSIM
The Structural SIMilarity (SSIM) index between two images, 图像结构相似性，是一种衡量两幅图像相似度的指标。用来判断图片压缩后的质量。由于这个图像指标很流行，百度搜索一大堆，而且MATLAB有函数可以调用，故不再详细介绍。接下来，只是给出我平时用SSIM评价红外图像质量的MATLAB代码：
How to use SSIM in MATLAB:
clc;
close all;
clear;
% A=imread('GABF1_dde_S4.jpg');
% ref=imread('S4.tif');
load('S1.mat');
load('S1_GABF.mat');
A= im2double(OUT);
ref = im2double(S);%S1-S4.mat
%ref = im2double(S11);%S11.mat
%调整图像大小有三种方法
%分别是近邻内插法、双线性内插法和双立方内插法
%使用的函数是imresize(src,size,method)
%A=imresize(A,[254,324],'nearest'); %近邻内插法,调整图像大小为254×324
%调整图像的维度
%ref=reshape(ref,254,324,3);
% ref1=ref(:,:,1);
%
% ref2=ref(:,:,2);
% ref3=ref(:,:,3);
% ref_a=zeros(254,324,3);%先申请变量空间;
% ref_a(:,:,1)=ref1;
% ref_a(:,:,2)=ref2;
% ref_a(:,:,3)=ref3;
% %ref=horzcat(ref1,ref2,ref3);
%
% A=imresize(A,[256,324],'bilinear'); %双线性内插法,调整图像大小为256×324
% %A=imresize(A,4,'bicubic'); %双立方内插法,放大4倍
% %ref ima
% ref_a=imresize(ref_a,[256,324],'bilinear'); %双线性内插法,调整图像大小为256×324
% ref=im2uint8(ref_a);
% ref=ref(:,:,3);
% A=A(:,:,3);
mval = ssim(A,ref)
上面的MATLAB代码，注释为以前的编写的，输入图像和参考图像以tif和bmp格式。没注释则为今天编写的，为mat格式，更加准确，一般经典的图像增强算法，输入图像都是mat格式。即把结果保存为mat格式，用save，直接在命令窗口输入help save查看保存为mat格式的具体用法。论文的算法，输入也是mat格式。
更新于9月3日，通用为了避免繁琐的操作，即评价多张图像的SSIM值，自然想到了函数和循环结构。故与上面的RMSC一致，还是直接在顶层调用SSIM，批量输出的方法：在循环里把多张图像用元胞向量的索引来访问输入图像，赋值给SSIM函数处理后，存储在矩阵向量中，遍历完所有图像后，计算平均值即可。
clc; close all; clear; %--------The reference image ref1=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/S1.png'); ref2=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/S2.png'); ref3=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/S3.png'); ref4=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/S4.png'); ref5=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/S11.png'); %-%--------The reference image for ssim %--------without the ROGABF % I1=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF1.png'); % I2=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF2.png'); % I3=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF3.png'); % I4=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF4.png'); % I5=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF11.png'); %--------without the ROGABF %--------without the wFilter I1=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter1.png'); I2=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter2.png'); I3=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter3.png'); I4=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter4.png'); I5=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_wFilter11.png'); %------without the wFilter %output five result of ssim at a time %the reference image ssim_out = zeros(1,5); ref = cell(1,5); ref{1,1}=ref1;ref{1,2}=ref2; ref{1,3}=ref3;ref{1,4}=ref4;ref{1,5}=ref5; %the reference image I = cell(1,5); I{1,1}=I1;I{1,2}=I2; I{1,3}=I3;I{1,4}=I4;I{1,5}=I5; for k = 1:5 tempA = double(I{k}); tempref = double(ref{k}); ssim_out(k) = ssim(tempA,tempref); end ssim_out mean(ssim_out) %output five result of RMSC at a time
重要的是想法：遇到多张图像或多个操作，能够想到函数加上循环结构，能够实现批量输出。如经常碰到的一次显示四种算法的实验结果，还有一次输出多张图像的仿真效果。
时隔半年后，再次更新，即原代码中把处理后图像与参考图像都转换为im2doubel型，在实验中发现SSIM值为0.2左右，与图像主观视觉相差极大，经过仔细看ssim的用法即输入数据类型，发现用im2uint8转换为uint8比较靠谱。
clc;
close all;
clear;
%--------The reference image
ref1=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/13raw01.png');
ref2=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/15seq4_nuc01.png');
ref3=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/04red_cmb01.png');
ref4=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/01orange_raw01.png');
ref5=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/19sempach_BG3_1.png');
%-%--------The reference image for ssim
%--------without the ROGABF
% I1=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF1.png');
% I2=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF2.png');
% I3=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF3.png');
% I4=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF4.png');
% I5=imread('E:\paper\MATLAB_CODE\GABF_DDE93_Ablation/ablation/SA_ROGABF11.png');
%--------without the ROGABF
%--------without the wFilter
%He
% I1=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/13raw01LD2.bmp');
I1=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/13raw01HE.png');
I2=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/15seq4_nuc01HE.png');
I3=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/04red_cmb01He.png');
I4=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/01orange_raw01He.png');
I5=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/19sempach_BG3_1He.png');
%He
%AHP_BC
% I1=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/IMA1LD.png');
% I1=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/13raw01AHP_BC1.bmp');
% I2=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/15seq4_nuc01AHP_BC.png');
% I3=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/04red_cmb01AHP_BC.png');
% I4=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/01orange_raw01AHP_BC.png');
% %324*256
% % %512*512
% I5=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/19sempach_BG3_1AHP_BC.png');
%AHP_BC
%proposed method
% I1=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/13raw01proposed1.bmp');
% I2=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/15seq4_nuc01HRM.png');
% I3=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/04red_cmb01HRLD.png');
% I4=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/01orange_raw01HRLD.png');
% %324*256
% %512*512
% I5=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/exp2/19sempach_BG3_1HRLD.png');
%proposed method modified
% I1=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/IMA1HRM.bmp');
% I2=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/IMA2HRM.bmp');
% I3=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/IMA3HRM.bmp');
% I4=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/IMA4HRM.bmp');
% I5=imread('E:/paper/MATLAB_CODE/MAT_fpga/pic/IMA5HRM.bmp');
%%proposed method modified
%output five result of ssim at a time
%the reference image
ssim_out = zeros(1,5);
ref = cell(1,5);
ref{1,1}=ref1;ref{1,2}=ref2;
ref{1,3}=ref3;ref{1,4}=ref4;ref{1,5}=ref5;
%the reference image
I = cell(1,5);
I{1,1}=I1;I{1,2}=I2;
I{1,3}=I3;I{1,4}=I4;I{1,5}=I5;
for k = 1:5
%------entropy计算-----------%
tempA = double(I{k});
tempA =uint8( (tempA));%the key of code1
ssim_out(k) =entropy(tempA);%entropy(I)
%------entropy计算-----------%
% %------SSIM计算-----------%
% %适用于图像I尺寸比参考图像REF大则调整图像I尺寸
% %图像为uint８型
% tempA = double(I{k});
% tempref = double(ref{k});
% [rm,rn]= size(tempref);
% tempA=imresize(tempA,[rm rn]);
% tempA = tempA(:,:,3);
% tempA =im2uint8( (tempA));%the key of code1
% tempref=im2uint8(255*mat2gray(tempref));%the key of code2
%
% % tempA = uint8(tempA);
% % % ssim_out(k) =entropy(tempA);%entropy(I)
% % tempA = uint8(tempA);
%
% ssim_out(k) = ssim(tempA,tempref);%ssim
% %------SSIM计算-----------%
% %-----PSNR计算-----------%
% %适用于图像I尺寸比参考图像REF大则调整图像I尺寸
% %图像为uint８型
% tempA = double(I{k});
% tempref = double(ref{k});
% [rm,rn]= size(tempref);
% tempA=imresize(tempA,[rm rn]);
% tempA = tempA(:,:,3);
% tempA =im2uint8( (tempA));%the key of code1
% tempref=im2uint8(255*mat2gray(tempref));%the key of code2
%
% % tempA = uint8(tempA);
% % % ssim_out(k) =entropy(tempA);%entropy(I)
% % tempA = uint8(tempA);
%
% ssim_out(k) = psnr(tempA,tempref);%ssim
% %------PSNR计算-----------%
end
ssim_out
mean(ssim_out)
%output five result of RMSC at a time
这次更新的代码包含了图像信息熵，SSIM和PSNR等场景的客观指标，在代码中用分割线标注，注意PSNR用im2uint8转换，其他用uint8转换，数据类型转换正确输出的结果为正常范围。
PSNR高于40dB说明图像质量极好（即非常接近原始图像），在30—40dB通常表示图像质量是好的（即失真可以察觉但可以接受），在20—30dB说明图像质量差；
最后，PSNR低于20dB图像不可接受。SSIM的正常范围为-1到1，越大表示与参考图像越相似。本文的图像信息熵都在1-10左右，暂未查找到有关信息熵正常范围的相关文献
三 FSIM
feature similarity index mersure(FSIM)，是SSIM的变种，该算法认为一张图片中的所有像素并非具有相同的重要性，比如物体边缘的像素点对于界定物体的结构肯定比其他背景区域的像素点更为重要；另外一种重要的评价指标VIF尽管在不同的子带上具有不同的权重，但是在具体的某一子带上参与计算的像素点均具有相同的权重；根据图像本身的特性，这样不加区分并不合适。因此改进的方向实际上重在如何区分这些重要点并给与合适的权重。在GitHub上可以搜到代码，要想进一步理解，可以去搜作者的论文。
在MATLAB软件里有函数可以调用，直接help FSIM, 即[FSIM, FSIMc] = FeatureSIM(img1, img2);
function [FSIM, FSIMc] = FSIM(imageRef, imageDis)
% ========================================================================
% FSIM Index with automatic downsampling, Version 1.0
% Copyright(c) 2010 Lin ZHANG, Lei Zhang, Xuanqin Mou and David Zhang
% All Rights Reserved.
%
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining 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.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for calculating the
% Feature SIMilarity (FSIM) index between two images.
%
% Please refer to the following paper
%
% Lin Zhang, Lei Zhang, Xuanqin Mou, and David Zhang,"FSIM: a feature similarity
% index for image qualtiy assessment", IEEE Transactions on Image Processing, vol. 20, no. 8, pp. 2378-2386, 2011.
%
%----------------------------------------------------------------------
%
%Input : (1) imageRef: the first image being compared
% (2) imageDis: the second image being compared
%
%Output: (1) FSIM: is the similarty score calculated using FSIM algorithm. FSIM
% only considers the luminance component of images. For colorful images,
% they will be converted to the grayscale at first.
% (2) FSIMc: is the similarity score calculated using FSIMc algorithm. FSIMc
% considers both the grayscale and the color information.
%Note: For grayscale images, the returned FSIM and FSIMc are the same.
%
%-----------------------------------------------------------------------
%
%Usage:
%Given 2 test images img1 and img2. For gray-scale images, their dynamic range should be 0-255.
%For colorful images, the dynamic range of each color channel should be 0-255.
%
%[FSIM, FSIMc] = FeatureSIM(img1, img2);
%-----------------------------------------------------------------------
[rows, cols] = size(imageRef(:,:,1));
I1 = ones(rows, cols);
I2 = ones(rows, cols);
Q1 = ones(rows, cols);
Q2 = ones(rows, cols);
if ndims(imageRef) == 3 %images are colorful
Y1 = 0.299 * double(imageRef(:,:,1)) + 0.587 * double(imageRef(:,:,2)) + 0.114 * double(imageRef(:,:,3));
Y2 = 0.299 * double(imageDis(:,:,1)) + 0.587 * double(imageDis(:,:,2)) + 0.114 * double(imageDis(:,:,3));
I1 = 0.596 * double(imageRef(:,:,1)) - 0.274 * double(imageRef(:,:,2)) - 0.322 * double(imageRef(:,:,3));
I2 = 0.596 * double(imageDis(:,:,1)) - 0.274 * double(imageDis(:,:,2)) - 0.322 * double(imageDis(:,:,3));
Q1 = 0.211 * double(imageRef(:,:,1)) - 0.523 * double(imageRef(:,:,2)) + 0.312 * double(imageRef(:,:,3));
Q2 = 0.211 * double(imageDis(:,:,1)) - 0.523 * double(imageDis(:,:,2)) + 0.312 * double(imageDis(:,:,3));
else %images are grayscale
Y1 = imageRef;
Y2 = imageDis;
end
Y1 = double(Y1);
Y2 = double(Y2);
%%%%%%%%%%%%%%%%%%%%%%%%%
% Downsample the image
%%%%%%%%%%%%%%%%%%%%%%%%%
minDimension = min(rows,cols);
F = max(1,round(minDimension / 256));
aveKernel = fspecial('average',F);
aveI1 = conv2(I1, aveKernel,'same');
aveI2 = conv2(I2, aveKernel,'same');
I1 = aveI1(1:F:rows,1:F:cols);
I2 = aveI2(1:F:rows,1:F:cols);
aveQ1 = conv2(Q1, aveKernel,'same');
aveQ2 = conv2(Q2, aveKernel,'same');
Q1 = aveQ1(1:F:rows,1:F:cols);
Q2 = aveQ2(1:F:rows,1:F:cols);
aveY1 = conv2(Y1, aveKernel,'same');
aveY2 = conv2(Y2, aveKernel,'same');
Y1 = aveY1(1:F:rows,1:F:cols);
Y2 = aveY2(1:F:rows,1:F:cols);
%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the phase congruency maps
%%%%%%%%%%%%%%%%%%%%%%%%%
PC1 = phasecong2(Y1);
PC2 = phasecong2(Y2);
%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the gradient map
%%%%%%%%%%%%%%%%%%%%%%%%%
dx = [3 0 -3; 10 0 -10; 3 0 -3]/16;
dy = [3 10 3; 0 0 0; -3 -10 -3]/16;
IxY1 = conv2(Y1, dx, 'same');
IyY1 = conv2(Y1, dy, 'same');
gradientMap1 = sqrt(IxY1.^2 + IyY1.^2);
IxY2 = conv2(Y2, dx, 'same');
IyY2 = conv2(Y2, dy, 'same');
gradientMap2 = sqrt(IxY2.^2 + IyY2.^2);
%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the FSIM
%%%%%%%%%%%%%%%%%%%%%%%%%
T1 = 0.85; %fixed
T2 = 160; %fixed
PCSimMatrix = (2 * PC1 .* PC2 + T1) ./ (PC1.^2 + PC2.^2 + T1);
gradientSimMatrix = (2*gradientMap1.*gradientMap2 + T2) ./(gradientMap1.^2 + gradientMap2.^2 + T2);
PCm = max(PC1, PC2);
SimMatrix = gradientSimMatrix .* PCSimMatrix .* PCm;
FSIM = sum(sum(SimMatrix)) / sum(sum(PCm));
%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the FSIMc
%%%%%%%%%%%%%%%%%%%%%%%%%
T3 = 200;
T4 = 200;
ISimMatrix = (2 * I1 .* I2 + T3) ./ (I1.^2 + I2.^2 + T3);
QSimMatrix = (2 * Q1 .* Q2 + T4) ./ (Q1.^2 + Q2.^2 + T4);
lambda = 0.03;
SimMatrixC = gradientSimMatrix .* PCSimMatrix .* real((ISimMatrix .* QSimMatrix) .^ lambda) .* PCm;
FSIMc = sum(sum(SimMatrixC)) / sum(sum(PCm));
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ResultPC]=phasecong2(im)
% ========================================================================
% Copyright (c) 1996-2009 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% The software is provided "as is", without warranty of any kind.
% References:
%
% Peter Kovesi, "Image Features From Phase Congruency". Videre: A
% Journal of Computer Vision Research. MIT Press. Volume 1, Number 3,
% Summer 1999 http://mitpress.mit.edu/e-journals/Videre/001/v13.html
nscale = 4; % Number of wavelet scales.
norient = 4; % Number of filter orientations.
minWaveLength = 6; % Wavelength of smallest scale filter.
mult = 2; % Scaling factor between successive filters.
sigmaOnf = 0.55; % Ratio of the standard deviation of the
% Gaussian describing the log Gabor filter's
% transfer function in the frequency domain
% to the filter center frequency.
dThetaOnSigma = 1.2; % Ratio of angular interval between filter orientations
% and the standard deviation of the angular Gaussian
% function used to construct filters in the
% freq. plane.
k = 2.0; % No of standard deviations of the noise
% energy beyond the mean at which we set the
% noise threshold point.
% below which phase congruency values get
% penalized.
epsilon = .0001; % Used to prevent division by zero.
thetaSigma = pi/norient/dThetaOnSigma; % Calculate the standard deviation of the
% angular Gaussian function used to
% construct filters in the freq. plane.
[rows,cols] = size(im);
imagefft = fft2(im); % Fourier transform of image
zero = zeros(rows,cols);
EO = cell(nscale, norient); % Array of convolution results.
estMeanE2n = [];
ifftFilterArray = cell(1,nscale); % Array of inverse FFTs of filters
% Pre-compute some stuff to speed up filter construction
% Set up X and Y matrices with ranges normalised to +/- 0.5
% The following code adjusts things appropriately for odd and even values
% of rows and columns.
if mod(cols,2)
xrange = [-(cols-1)/2:(cols-1)/2]/(cols-1);
else
xrange = [-cols/2:(cols/2-1)]/cols;
end
if mod(rows,2)
yrange = [-(rows-1)/2:(rows-1)/2]/(rows-1);
else
yrange = [-rows/2:(rows/2-1)]/rows;
end
[x,y] = meshgrid(xrange, yrange);
radius = sqrt(x.^2 + y.^2); % Matrix values contain *normalised* radius from centre.
theta = atan2(-y,x); % Matrix values contain polar angle.
% (note -ve y is used to give +ve
% anti-clockwise angles)
radius = ifftshift(radius); % Quadrant shift radius and theta so that filters
theta = ifftshift(theta); % are constructed with 0 frequency at the corners.
radius(1,1) = 1; % Get rid of the 0 radius value at the 0
% frequency point (now at top-left corner)
% so that taking the log of the radius will
% not cause trouble.
sintheta = sin(theta);
costheta = cos(theta);
clear x; clear y; clear theta; % save a little memory
% Filters are constructed in terms of two components.
% 1) The radial component, which controls the frequency band that the filter
% responds to
% 2) The angular component, which controls the orientation that the filter
% responds to.
% The two components are multiplied together to construct the overall filter.
% Construct the radial filter components...
% First construct a low-pass filter that is as large as possible, yet falls
% away to zero at the boundaries. All log Gabor filters are multiplied by
% this to ensure no extra frequencies at the 'corners' of the FFT are
% incorporated as this seems to upset the normalisation process when
% calculating phase congrunecy.
lp = lowpassfilter([rows,cols],.45,15); % Radius .45, 'sharpness' 15
logGabor = cell(1,nscale);
for s = 1:nscale
wavelength = minWaveLength*mult^(s-1);
fo = 1.0/wavelength; % Centre frequency of filter.
logGabor{s} = exp((-(log(radius/fo)).^2) / (2 * log(sigmaOnf)^2));
logGabor{s} = logGabor{s}.*lp; % Apply low-pass filter
logGabor{s}(1,1) = 0; % Set the value at the 0 frequency point of the filter
% back to zero (undo the radius fudge).
end
% Then construct the angular filter components...
spread = cell(1,norient);
for o = 1:norient
angl = (o-1)*pi/norient; % Filter angle.
% For each point in the filter matrix calculate the angular distance from
% the specified filter orientation. To overcome the angular wrap-around
% problem sine difference and cosine difference values are first computed
% and then the atan2 function is used to determine angular distance.
ds = sintheta * cos(angl) - costheta * sin(angl); % Difference in sine.
dc = costheta * cos(angl) + sintheta * sin(angl); % Difference in cosine.
dtheta = abs(atan2(ds,dc)); % Absolute angular distance.
spread{o} = exp((-dtheta.^2) / (2 * thetaSigma^2)); % Calculate the
% angular filter component.
end
% The main loop...
EnergyAll(rows,cols) = 0;
AnAll(rows,cols) = 0;
for o = 1:norient % For each orientation.
sumE_ThisOrient = zero; % Initialize accumulator matrices.
sumO_ThisOrient = zero;
sumAn_ThisOrient = zero;
Energy = zero;
for s = 1:nscale, % For each scale.
filter = logGabor{s} .* spread{o}; % Multiply radial and angular
% components to get the filter.
ifftFilt = real(ifft2(filter))*sqrt(rows*cols); % Note rescaling to match power
ifftFilterArray{s} = ifftFilt; % record ifft2 of filter
% Convolve image with even and odd filters returning the result in EO
EO{s,o} = ifft2(imagefft .* filter);
An = abs(EO{s,o}); % Amplitude of even & odd filter response.
sumAn_ThisOrient = sumAn_ThisOrient + An; % Sum of amplitude responses.
sumE_ThisOrient = sumE_ThisOrient + real(EO{s,o}); % Sum of even filter convolution results.
sumO_ThisOrient = sumO_ThisOrient + imag(EO{s,o}); % Sum of odd filter convolution results.
if s==1 % Record mean squared filter value at smallest
EM_n = sum(sum(filter.^2)); % scale. This is used for noise estimation.
maxAn = An; % Record the maximum An over all scales.
else
maxAn = max(maxAn, An);
end
end % ... and process the next scale
% Get weighted mean filter response vector, this gives the weighted mean
% phase angle.
XEnergy = sqrt(sumE_ThisOrient.^2 + sumO_ThisOrient.^2) + epsilon;
MeanE = sumE_ThisOrient ./ XEnergy;
MeanO = sumO_ThisOrient ./ XEnergy;
% Now calculate An(cos(phase_deviation) - | sin(phase_deviation)) | by
% using dot and cross products between the weighted mean filter response
% vector and the individual filter response vectors at each scale. This
% quantity is phase congruency multiplied by An, which we call energy.
for s = 1:nscale,
E = real(EO{s,o}); O = imag(EO{s,o}); % Extract even and odd
% convolution results.
Energy = Energy + E.*MeanE + O.*MeanO - abs(E.*MeanO - O.*MeanE);
end
% Compensate for noise
% We estimate the noise power from the energy squared response at the
% smallest scale. If the noise is Gaussian the energy squared will have a
% Chi-squared 2DOF pdf. We calculate the median energy squared response
% as this is a robust statistic. From this we estimate the mean.
% The estimate of noise power is obtained by dividing the mean squared
% energy value by the mean squared filter value
medianE2n = median(reshape(abs(EO{1,o}).^2,1,rows*cols));
meanE2n = -medianE2n/log(0.5);
estMeanE2n(o) = meanE2n;
noisePower = meanE2n/EM_n; % Estimate of noise power.
% Now estimate the total energy^2 due to noise
% Estimate for sum(An^2) + sum(Ai.*Aj.*(cphi.*cphj + sphi.*sphj))
EstSumAn2 = zero;
for s = 1:nscale
EstSumAn2 = EstSumAn2 + ifftFilterArray{s}.^2;
end
EstSumAiAj = zero;
for si = 1:(nscale-1)
for sj = (si+1):nscale
EstSumAiAj = EstSumAiAj + ifftFilterArray{si}.*ifftFilterArray{sj};
end
end
sumEstSumAn2 = sum(sum(EstSumAn2));
sumEstSumAiAj = sum(sum(EstSumAiAj));
EstNoiseEnergy2 = 2*noisePower*sumEstSumAn2 + 4*noisePower*sumEstSumAiAj;
tau = sqrt(EstNoiseEnergy2/2); % Rayleigh parameter
EstNoiseEnergy = tau*sqrt(pi/2); % Expected value of noise energy
EstNoiseEnergySigma = sqrt( (2-pi/2)*tau^2 );
T = EstNoiseEnergy + k*EstNoiseEnergySigma; % Noise threshold
% The estimated noise effect calculated above is only valid for the PC_1 measure.
% The PC_2 measure does not lend itself readily to the same analysis. However
% empirically it seems that the noise effect is overestimated roughly by a factor
% of 1.7 for the filter parameters used here.
T = T/1.7; % Empirical rescaling of the estimated noise effect to
% suit the PC_2 phase congruency measure
Energy = max(Energy - T, zero); % Apply noise threshold
EnergyAll = EnergyAll + Energy;
AnAll = AnAll + sumAn_ThisOrient;
end % For each orientation
ResultPC = EnergyAll ./ AnAll;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOWPASSFILTER - Constructs a low-pass butterworth filter.
%
% usage: f = lowpassfilter(sze, cutoff, n)
%
% where: sze is a two element vector specifying the size of filter
% to construct [rows cols].
% cutoff is the cutoff frequency of the filter 0 - 0.5
% n is the order of the filter, the higher n is the sharper
% the transition is. (n must be an integer >= 1).
% Note that n is doubled so that it is always an even integer.
%
% 1
% f = --------------------
% 2n
% 1.0 + (w/cutoff)
%
% The frequency origin of the returned filter is at the corners.
%
% See also: HIGHPASSFILTER, HIGHBOOSTFILTER, BANDPASSFILTER
%
% Copyright (c) 1999 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% October 1999
% August 2005 - Fixed up frequency ranges for odd and even sized filters
% (previous code was a bit approximate)
function f = lowpassfilter(sze, cutoff, n)
if cutoff < 0 || cutoff > 0.5
error('cutoff frequency must be between 0 and 0.5');
end
if rem(n,1) ~= 0 || n < 1
error('n must be an integer >= 1');
end
if length(sze) == 1
rows = sze; cols = sze;
else
rows = sze(1); cols = sze(2);
end
% Set up X and Y matrices with ranges normalised to +/- 0.5
% The following code adjusts things appropriately for odd and even values
% of rows and columns.
if mod(cols,2)
xrange = [-(cols-1)/2:(cols-1)/2]/(cols-1);
else
xrange = [-cols/2:(cols/2-1)]/cols;
end
if mod(rows,2)
yrange = [-(rows-1)/2:(rows-1)/2]/(rows-1);
else
yrange = [-rows/2:(rows/2-1)]/rows;
end
[x,y] = meshgrid(xrange, yrange);
radius = sqrt(x.^2 + y.^2); % A matrix with every pixel = radius relative to centre.
f = ifftshift( 1 ./ (1.0 + (radius ./ cutoff).^(2*n)) ); % The filter
return;
上面的MATLAB代码为我在网上搜索的FSIM代码，还未详细研究和调用，故以后再进一步介绍。
四 MOS
the mean opinion score(MOS)，一种主观的评价方式，即对图像的处理效果打分，我的理解为观察者对算法的处理结果从对比度，真实性，细节清晰度等三方面进行打分，如每一项满分5分，值越大说明处理结果越好，具体的介绍见下面的文献[ref]：
【ref】Realtime infrared image detail enhancement based on fast guided image filter and plateau equalization.
五 信息熵（entropy）
图像中平均信息量的多少。这是一种简单的评价方式，MATLAB代码在网上和matlab软件里都有函数可以调用。
The matlab code of entropy
%求图像熵值
%传入一张彩色图片的矩阵
%输出图片的图像熵值
I_gray=imread('S11_ROG1.tif');
%I_gray = rgb2gray(I);
[ROW,COL] = size(I_gray);
%%
%新建一个size =256的矩阵，用于统计256个灰度值的出现次数
temp = zeros(256);
for i= 1:ROW
for j = 1:COL
%统计当前灰度出现的次数
temp(I_gray(i,j)+1)= temp(I_gray(i,j)+1)+1;
end
end
%%
res = 0.0 ;
for i = 1:256
%计算当前灰度值出现的概率
temp(i) = temp(i)/(ROW*COL);
%如果当前灰度值出现的次数不为0
if temp(i)~=0.0
res = res - temp(i) * (log(temp(i)) / log(2.0));
end
end
disp(res);
当然图像的评价指标还有psnr, mse等都是常见的，由于简单，故不会作为论文中的客观评价指标。本文选取的为RMSC, SSIM, MOS三种。还有算法运行时间，这个与MATLAB软件，个人电脑的内存，CPU型号，系统等有关，故作为评价指标并不可靠，谨记！通常用算法复杂度衡量更准确。算法运行时间，在matlab中用tic与 toc函数分别加在算法的开头与末尾即可。
六 常见图像处理指标（图像的均值和标准差）
这些基本的代码是得掌握，直接贴出matlab代码如下：
1. 求图像的均值
close all;
clear;
clc;
i=imread('d:/lena.jpg'); %载入真彩色图像路径
i=rgb2gray(i); %转换为灰度图
i=double(i); %将uint8型转换为double型，否则不能计算统计量
[m,n]=size(i);
s=0;
for x=1:m
for y=1:n
s=s+i(x,y); %求像素值总和 s ， i(x,y)表示位于某个坐标下的像素值
end
end
%所有像素均值
a1=mean(mean(i)); %第一种方法：先计算列向量均值，再求总均值。
a2=mean2(i); %第二种方法：用函数mean2求总均值
a3=s/(m*n); %第三种方法：按公式计算，像素值总和除以像素个数。
a4=sum(sum(i))/(m*n); %第四种方法：也是按公式计算，但是用sum来求像素值总和。
2.求图像的标准差
close all
clear
clc;
i=imread('d:/lena.jpg'); %载入真彩色图像
i=rgb2gray(i); %转换为灰度图
i=double(i); %将uint8型转换为double型，否则不能计算统计量
avg=mean2(i); %求图像均值
[m,n]=size(i);
s=0;
for x=1:m
for y=1:n
s=s+(i(x,y)-avg)^2; %求得所有像素与均值的平方和。
end
end
%求图像的方差
a1=var(i(:)); %第一种方法：利用函数var求得。
a2=s/(m*n-1); %第二种方法：利用方差公式求得
a3=(std2(i))^2; %第三种方法：利用std2求得标准差，再平方即为方差。
六 MATLAB调试（找错与改正）
错误使用 conv2
不支持 N 维数组。
出错 filter2 (line 59)第二级
y = conv2(hcol, hrow, x, shape);
出错 ssim (line 188)第一级
mu1 = filter2(window, img1, 'valid');
出错 SSIM_CAL (line 24)顶层
mval = ssim(A,ref)
%%上面为报错信息
1.找错:
一是看懂MATLAB提示，最底下为最顶层函数1，往上为第一级函数1.1，第二级函数1.1.1。在上面的例子，SSIM_CAL 为顶层函数，在其24行的ssim函数（第一级）出错。在ssim函数的188行的filter2函数（第二级）出错。再继续深入到filter2函数中的59行出错，即conv2函数出错；就是问题源头所在、已找到错误。
2.改正（解决问题）
设置断点，即在问题源头出设置断点，运行matlab程序，然后再查看断点处函数的值或变量值，找到出错原因，并改正。比如，本例子中则在filter2函数中的59行（问题源头）设置断点，运行后查看这一行代码conv2的变量值，和conv2函数的用法和注意事项，即conv2函数仅支持2维数组卷积，但运行到这个断点，发现conv2函数中的x变量，即输入图像为m*n*3，即多维数组，故报错。通过断点找到原因后，解决问题十分容易，即将输入图像搞成二维即可。
这就是调试matlab程序的方法，通过报错定位问题，再通过设置断点找出问题，再根据问题源头解决即可！
3.函数+循环结构= 飞一般的感觉
遇到输入多张图像或多个数据，马上想到函数加循环，别手动慢慢跑仿真了，太耗时。只要是多个的循环操作都能用函数加循环解决。刚开始编写批量输出的顶层函数，有点费劲，当你打通了，解决这种多个操作和数据问题的思路就很不一样，那么后面就大大节省时间，因为思路相同，不同的问题套进去就行。
想法是：
安逸和享乐绝不是生活的追求！
-----------爱因斯坦
朝最大抵抗力走！
-----------朱光潜
这是今天做红外图像细节增强算法的收获！
