基于SSIM （结构相似性）的图像比较

function [mssim, ssim_map,siga_sq,sigb_sq] = ssim2(ima, imb)  
% ========================================================================  
%ssim的算法主要参考如下论文：  
%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.  
%  首先对图像加窗处理，w=fspecial('gaussian', 11, 1.5);  
%                 (2*ua*ub+C1)*(2*sigmaa*sigmab+C2)  
%   SSIM(A,B)=————————————————————————  
%              (ua*ua+ub*ub+C1)(sigmaa*sigmaa+sigmab*sigmab+C2)  
%     C1=（K1*L）;  
%     C2=(K2*L);   K1=0.01，K2=0.03  
%     L为灰度级数，L=255  
%-------------------------------------------------------------------  
%     ima - 比较图像A  
%     imb - 比较图像B  
%  
% ssim_map - 各加窗后得到的SSIM（A,B|w）组成的映射矩阵  
%    mssim - 对加窗得到的SSIM（A,B|w）求平均，即最终的SSIM（A,B）  
%  siga_sq - 图像A各窗口内灰度值的方差  
%  sigb_sq - 图像B各窗口内灰度值的方差  
%-------------------------------------------------------------------  
%  Cool_ben  
%========================================================================  

w = fspecial('gaussian', 11, 1.5);  %window 加窗  
K(1) = 0.01;                      
K(2) = 0.03;                      
L = 255;       
ima = double(ima);  
imb = double(imb);  
tic
C1 = (K(1)*L)^2;  
C2 = (K(2)*L)^2;  
w = w/sum(sum(w));  

ua   = filter2(w, ima, 'valid');%对窗口内并没有进行平均处理，而是与高斯卷积，  
ub   = filter2(w, imb, 'valid'); % 类似加权平均  
ua_sq = ua.*ua;  
ub_sq = ub.*ub;  
ua_ub = ua.*ub;  
siga_sq = filter2(w, ima.*ima, 'valid') - ua_sq;  
sigb_sq = filter2(w, imb.*imb, 'valid') - ub_sq;  
sigab = filter2(w, ima.*imb, 'valid') - ua_ub;  

ssim_map = ((2*ua_ub + C1).*(2*sigab + C2))./((ua_sq + ub_sq + C1).*(siga_sq + sigb_sq + C2));  
toc

mssim = mean2(ssim_map);  

return  

测试代码如下：

clearvars;
close all;
clc;
img = imread('D:/picture/lena.jpg');
gray = rgb2gray(img);
op = fspecial('average',3);
gray_ave = imfilter(gray, op,'replicate');
[mssim2, ssim_map2,siga_sq2,sigb_sq2] = ssim2(gray, gray_ave);  
[minlocy, minlocx] = find(ssim_map2 == min(ssim_map2(:)));
[maxlocy, maxlocx] = find(ssim_map2 == max(ssim_map2(:)));
figure, imshow(ssim_map2); 
hold on;
plot(minlocx, minlocy, 'rs');
plot(maxlocx, maxlocy, 'go');
hold off;
figure, imshow(gray); title('original');
hold on;
plot(minlocx, minlocy, 'rs');
plot(maxlocx, maxlocy, 'go');
hold off;
figure, imshow(gray_ave); title('average');
hold on;
plot(minlocx, minlocy, 'rs');
plot(maxlocx, maxlocy, 'go');
hold off;