计算matlab中图像的PSNR和SSIM

Dear 丶 2022-05-08 11:50 355阅读 0赞

网上找了很多关于PSNR和SSIM的计算，很多结果算出来都不一样，公式都是普遍的，如下：

![Center][]

现在总结下造成结果差异的原因。

PSNR的差异：

1.灰度图像：灰度图像比较好计算，只有一个灰度值。

2.彩色图像：

**(a)**可以将分别计算R,G,B三个通道总和，最后MSE直接在原公式上多除以3就行（[opencv官方代][opencv][码][opencv]是这么做的，与[matlab直接计算结果][matlab]是一样的）。

**(b)**将R,G,B格式转换为YCbCr，只计算Y分量（亮度分量），结果会比直接计算要高几个dB。

贴代码，这里是将图片格式转成YCbCr（只计算Y分量）：

function [PSNR, MSE] = psnr(X, Y) 
       %%%%%%%%%%%%%%%%%%%%%%%%%%% 
       % 
       % 计算峰值信噪比PSNR 
       % 将RGB转成YCbCr格式进行计算 
       % 如果直接计算会比转后计算值要小 2dB左右（当然是个别测试） 
       % 
       %%%%%%%%%%%%%%%%%%%%%%%%%%% 
       if size(X, 3)~= 1 %判断图像时不是彩色图，如果是，结果为 3，否则为 1 
       org=rgb2ycbcr(X); 
       test=rgb2ycbcr(Y); 
       Y1=org(:,:, 1); 
       Y2=test(:,:, 1); 
       Y1= double(Y1); %计算平方时候需要转成 double类型，否则uchar类型会丢失数据 
       Y2= double(Y2); 
       else %灰度图像，不用转换 
       Y1= double(X); 
       Y2= double(Y); 
       end 
       
       if nargin< 2 
       D = Y1; 
       else 
       if any(size(Y1)~=size(Y2)) 
       error( 'The input size is not equal to each other!'); 
       end 
       D = Y1 - Y2; 
       end 
       MSE = sum(D(:).*D(:)) / numel(Y1); 
       PSNR = 10* log10( 255^ 2 / MSE);

控制台输入下面三条语句：

>> X= imread( 'C:\Users\Administrator\Desktop\noise_image.jpg'); 
       >> Y= imread( 'C:\Users\Administrator\Desktop\actruel_image.jpg'); 
       >> psnr(X, Y)

SSIM的差异：同上，如果直接不转换成YCbCr格式，结果会偏高很多（ matlab中，SSIM提出者【1】， [代码][matlab] ）。opencv里面是分别计算了R,G,B三个分量的SSIM值（ [官方代码][opencv] ）。最后我将3个值取了个平均（这个值比matlab里面低很多）。

以下代码主要是参考原作者修改的，源代码是直接没有进行格式转换，直接RGB格式，下面我是将他转换成YCbCr计算图片的SSIM

function [mssim, ssim_map] = ssim(img1, img2, K, window, L) 
       
       %======================================================================== 
       %SSIM Index, Version 1.0 
       %Copyright(c) 2003 Zhou Wang 
       %All Rights Reserved. 
       % 
       %The author is with Howard Hughes Medical Institute, and Laboratory 
       % for Computational Vision at Center for Neural Science and Courant 
       %Institute of Mathematical Sciences, New York University. 
       % 
       %---------------------------------------------------------------------- 
       %Permission to use, copy, or modify 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 ap pearon 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 
       %Structural SIMilarity (SSIM) index between two images. Please refer 
       %to the following paper: 
       % 
       %Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image 
       %quality assessment: From error visibility to structural similarity" 
       %IEEE Transactios on Image Processing, vol. 13, no. 4, pp .600 -612, 
       %Apr. 2004. 
       % 
       %Kindly report any suggestions or corrections to zhouwang@ieee.org 
       % 
       %---------------------------------------------------------------------- 
       % 
       %Input : ( 1) img1: the first image being compared 
       % ( 2) img2: the second image being compared 
       % ( 3) K: constants in the SSIM index formula (see the above 
       % reference). defualt value: K = [ 0.01 0.03] 
       % ( 4) window: local window for statistics (see the above 
       % reference). default widnow is Gaussian given by 
       % window = fspecial( 'gaussian', 11, 1.5); 
       % ( 5) L: dynamic range of the images. default: L = 255 
       % 
       %Output: ( 1) mssim: the mean SSIM index value between 2 images. 
       % If one of the images being compared is regarded as 
       % perfect quality, then mssim can be considered as the 
       % quality measure of the other image. 
       % If img1 = img2, then mssim = 1. 
       % ( 2) ssim_map: the SSIM index map of the test image. The map 
       % has a smaller size than the input images. The actual size: 
       % size(img1) - size(window) + 1. 
       % 
       %Default Usage: 
       % Given 2 test images img1 and img2, whose dynamic range is 0 -255 
       % 
       % [mssim ssim_map] = ssim_index(img1, img2); 
       % 
       %Advanced Usage: 
       % User defined parameters. For example 
       % 
       % K = [ 0.05 0.05]; 
       % window = ones( 8); 
       % L = 100; 
       % [mssim ssim_map] = ssim_index(img1, img2, K, window, L); 
       % 
       %See the results: 
       % 
       % mssim %Gives the mssim value 
       % imshow(max( 0, ssim_map).^ 4) %Shows the SSIM index map 
       % 
       %======================================================================== 
       
       
       if (nargin < 2 | nargin > 5) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       
       if (size(img1) ~= size(img2)) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       
       [M N] = size(img1); 
       
       if (nargin == 2) 
       if ((M < 11) | (N < 11)) % 图像大小过小，则没有意义。 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return 
       end 
       window = fspecial( 'gaussian', 11, 1.5); % 参数一个标准偏差 1.5， 11* 11的高斯低通滤波。 
       K( 1) = 0.01; % default settings 
       K (2) = 0.03; 
       L = 255; 
       end 
       
       if (nargin == 3) 
       if ((M < 11) | (N < 11)) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return 
       end 
       window = fspecial( 'gaussian', 11, 1.5); 
       L = 255; 
       if (length(K) == 2) 
       if (K( 1) < 0 | K( 2) < 0) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       else 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       end 
       
       if (nargin == 4) 
       [H W] = size(window); 
       if ((H*W) < 4 | (H > M) | (W > N)) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return 
       end 
       L = 255; 
       if (length(K) == 2) 
       if (K( 1) < 0 | K( 2) < 0) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       else 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       end 
       
       if (nargin == 5) 
       [H W] = size(window); 
       if ((H*W) < 4 | (H > M) | (W > N)) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return 
       end 
       if (length(K) == 2) 
       if (K( 1) < 0 | K( 2) < 0) 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       else 
       ssim_index = -Inf; 
       ssim_map = -Inf; 
       return; 
       end 
       end 
       
       if size(img1,3)~= 1 %判断图像时不是彩色图，如果是，结果为 3，否则为 1 
       org=rgb2ycbcr(img1); 
       test=rgb2ycbcr(img2); 
       y1=org(:,:, 1); 
       y2=test(:,:, 1); 
       y1= double(y1); 
       y2= double(y2); 
       else 
       y1= double(img1); 
       y2= double(img2); 
       end 
       img1 = double(y1); 
       img2 = double(y2); 
       % 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; % 计算C1参数，给亮度L（x，y）用。 C1= 6.502500 
       C2 = (K( 2)*L)^ 2; % 计算C2参数，给对比度C（x，y）用。 C2= 58.522500 
       window = window/sum(sum(window)); %滤波器归一化操作。 
       
       
       mu1 = filter2(window, img1, 'valid'); % 对图像进行滤波因子加权 valid改成same结果会低一丢丢 
       mu2 = filter2(window, img2, 'valid'); % 对图像进行滤波因子加权 
       
       mu1_sq = mu1.*mu1; % 计算出Ux平方值。 
       mu2_sq = mu2.*mu2; % 计算出Uy平方值。 
       mu1_mu2 = mu1.*mu2; % 计算Ux*Uy值。 
       
       sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq; % 计算sigmax （标准差） 
       sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq; % 计算sigmay （标准差） 
       sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2; % 计算sigmaxy（标准差） 
       
       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 
       mssim = mean2(ssim_map); 
       
       return

控制台输入以下代码：

>> img1= imread( 'C:\Users\Administrator\Desktop\noise_image.jpg'); 
       >> img2= imread( 'C:\Users\Administrator\Desktop\actruel_image.jpg'); 
       >> ssim(img1,img2)

最后说一句，不管是结果如何，只要对比实验用的同一种评价代码工具，无所谓结果和原论文一不一样，问题是很多论文实验都搞不出来滴

参考文献

【1】Wang Z, Bovik A C, Sheikh H R, et al. Image quality assessment: from error visibility to structural similarity\[J\]. IEEE Transactions on Image Processing, 2004, 13(4):600-612.

[Center]: /images/20220508/778e13f5a54b40249253e80f3bc38beb.png
[opencv]: http://www.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/highgui/video-input-psnr-ssim/video-input-psnr-ssim.html?highlight=ssim
[matlab]: http://www.cns.nyu.edu/~lcv/ssim/ssim_index.m