PSNR计算

峰值信噪比（Peak Signal to Noise Ratio,PSNR）
1、是使用广泛的评价图像质量的标准之一；
2、应用于图像压缩，PSNR越大，表示图像的失真度就越小；
3、如果重构后的图像和原图像完全相同，则PSNR=Inf（无穷大），原图像为全白，重构图像为全黑，则PSNR=0；
4、通常对于有损压缩来说，PSNR值通常介于30~50dB之间；
5、对于使用相同的原图像和相同的压缩比的不同算法，PSNR越大，通常会说明该算法对图像的保真效果越好。

PSNR=10log10f2max∑Mx=1∑Ny=1(f(x,y)−f^(x,y))2

f(x,y)：表示原始图像；
f^(x,y):重构图像
M,N：图像的水平方向和垂直方向的宽度和高度
fmax:原始图像和重构图像中可能的最大像素值，对于8位灰度，该值为255

%彩色图像的PSNR计算 
% demo for PSNR function 
% Copyright (c) 2016,  
% Sunlinju 
% All Rights Reserved 
%1.the two image have the same size; 
%2.if  psnr=1000,it means the two image are same. 
clc; 
close all; 
clear all; 
imgOriginal=double(imread('DemCap_240G_On.bmp'));%must be double ,or not sumSqrError always equals 255 
imgProcessed=double(imread('DemCap_240G_Off.bmp')); 
height = size(imgOriginal,1); 
width = size(imgOriginal,2); 
errorR=0; 
errorG=0; 
errorB=0; 
sumSqrError=0; 
mse=0; 
psnr=0; 
for i=1:height 
    for j=1:width 
        errorR=imgProcessed(i,j,1)-imgOriginal(i,j,1); 
        errorG=imgProcessed(i,j,2)-imgOriginal(i,j,2); 
        errorB=imgProcessed(i,j,3)-imgOriginal(i,j,3); 
        sumSqrError=sumSqrError+errorR*errorR+errorG*errorG+errorB*errorB; 
    end 
end 
if sumSqrError~=0 
    mse = sumSqrError/(height*width*3); 
    psnr=10*log10(255*255/mse) 
else 
    psnr=1000  %means psnr=Inf 
end 

%波峰波谷检测 
%%  Peak  Test1 
% a = 20; 
% x = -30:1:30; 
% y = 1./sqrt(a^2+(a+x).^2)+1./sqrt(a^2+(a-x).^2); 
% plot(x,y); 
% grid on; 
% yy = diff(y); 
% yy(yy<0) = -1; 
% yy(yy>0) = 1; 
% yyy = diff(yy); 
%  
% mv = yyy(yyy~=0); 
% id = find(yyy~=0); 
% x0 = x(id); 
% y0 = y(id); 
% hold on; 
% for k = 1:length(id); 
%     plot(x0(k),y0(k),’r.’); 
%     if mv(k)<0 
%         s = ‘峰:’; 
%     else 
%         s = ‘谷:’; 
%     end; 
%     s = [s num2str(x0(k)) ‘,’ num2str(y0(k))]; 
%     text(x0(k),y0(k),s); 
% end; 
% hold off; 
%% Peak test2 
% data=[105.03    99.18    84.965    72.445    68.994    77.265    91.052    100.61    98.215    86.363    74.439    71.625    80.061    92.18    97.823    91.483    80.241    73.616    78.547    89.084    94.924    89.689    79.898    75.485    81.544    89.485    90.578    83.712    77.401    80.18    86.904    88.721    83.468    78.971    81.983    86.25    85.224    80.901    80.808    84.488] 
% IndMin=find(diff(sign(diff(data)))>0)+1; 
% IndMax=find(diff(sign(diff(data)))<0)+1; 
% plot(1:length(data),data); 
% hold on;grid on; 
% plot(IndMin,data(IndMin),’r^’) 
% plot(IndMax,data(IndMax),’k*’) 
% legend(‘曲线’,’波谷点’,’波峰点’) 
%% Peak test3 
% data=[ 2  12  4  6  9  4  3  1  19  7]; 
% plot(1:length(data),data); 
% hold on;grid on; 
% pks=findpeaks(data) %pks=  12   9  19 
%% 将Nx1宫格内数据归一化为其均值 
image_in=double(imread(‘DemCap_016R_Off.bmp’));%读取图像 
height = size(image_in,1); 
width = size(image_in,2); 
N=20; 
Template=zeros(1,N);

for i=1:N:height-N
for j=1:width
TemplateAvg=mean(mean(image_in(i:i+N-1,j,1)));
image_in(i:i+N-1,j,1)=TemplateAvg;
end
end
imshow(uint8(image_in));