圖像品質客觀評價名額
在做紅外圖像細節增強算法研究時，很重要一點就是要對經過算法處理的圖像結果進行評價，分成兩種評價方法。一種是視覺效果評價：主觀的人眼觀察，主要是通過觀察者能否看到更多圖像細節，給人的感覺是否真實自然，與其他算法處理的圖像的亮度和對比度的比較。在論文中的呈現方式，就是把各個算法的處理結果和原圖一起呈現給觀察者，再對圖像某一具體區域進行詳細分析，如圖像中的文字，黑暗區域，欄杆樹木和建築物等細節資訊進行比較，看那種算法的結果更清晰，暴露出更多的細節。當然咯，主觀評價方式往往取決于觀察者，并不是一種客觀可靠的方式，往往需要幾種公認的圖像品質客觀評價名額，對紅外圖像細節增強算法的處理結果進行比較。通過閱讀大量的SCI文獻，如IEEE Transactions on Image Processing，Optic Engineering, Applied Optics, Optik, Infrared Physics & Technology等多種權威期刊的圖像處理論文精讀，本文選取了四種經典的客觀評價名額：RMSC, SSIM, FSIM, MOS
一 RMSC
The Root-Mean-Square Contrast (RMSC) index，即均方根對比度，用來衡量圖像的對比度，其值越大對比度越高。其中， 為圖像總像素的均值。I為圖像某點的像素。M和N分别為圖像的高和寬。
The MATLAB code of RMSC is following:
close all
clear
clc;
load(‘S11.mat’);
I = double(S11);
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=sqrt(s/(m*n)); %利用方差公式求得. sqrt傳回數組 X 的每個元素的平方根 MATLAB Code
二 SSIM
The Structural SIMilarity (SSIM) index between two images, 圖像結構相似性，是一種衡量兩幅圖像相似度的名額。用來判斷圖檔壓縮後的品質。由于這個圖像名額很流行，百度搜尋一大堆，而且MATLAB有函數可以調用，故不再詳細介紹。接下來，隻是給出我平時用SSIM評價紅外圖像品質的MATLAB代碼：
How to use RMSC 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) The use of SSIM Code
上面的MATLAB代碼，注釋為以前的編寫的，輸入圖像和參考圖像以tif和bmp格式。沒注釋則為今天編寫的，為mat格式，更加準确，一般經典的圖像增強算法，輸入圖像都是mat格式。即把結果儲存為mat格式，用save，直接在指令視窗輸入help save檢視儲存為mat格式的具體用法。論文的算法，輸入也是mat格式。
三 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; FSIM Code
上面的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); Entropy Code
當然圖像的評價名額還有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程式的方法，通過報錯定位問題，再通過設定斷點找出問題，再根據問題源頭解決即可！
這是今天做紅外圖像細節增強算法的收獲！