ivan@119: function [mssim, ssim_map] = SMALL_ssim_index(img1, img2, K, window, L)
ivan@128: %% SSIM Index - Removed dependence on Image Processing Toolbox 
ivan@119: %========================================================================
ivan@119: %SSIM Index, Version 1.0
ivan@119: %Copyright(c) 2003 Zhou Wang
ivan@119: %All Rights Reserved.
ivan@119: %
ivan@119: %The author is with Howard Hughes Medical Institute, and Laboratory
ivan@119: %for Computational Vision at Center for Neural Science and Courant
ivan@119: %Institute of Mathematical Sciences, New York University.
ivan@119: %
ivan@119: %----------------------------------------------------------------------
ivan@119: %Permission to use, copy, or modify this software and its documentation
ivan@119: %for educational and research purposes only and without fee is hereby
ivan@119: %granted, provided that this copyright notice and the original authors'
ivan@119: %names appear on all copies and supporting documentation. This program
ivan@119: %shall not be used, rewritten, or adapted as the basis of a commercial
ivan@119: %software or hardware product without first obtaining permission of the
ivan@119: %authors. The authors make no representations about the suitability of
ivan@119: %this software for any purpose. It is provided "as is" without express
ivan@119: %or implied warranty.
ivan@119: %----------------------------------------------------------------------
ivan@119: %
ivan@119: %This is an implementation of the algorithm for calculating the
ivan@119: %Structural SIMilarity (SSIM) index between two images. Please refer
ivan@119: %to the following paper:
ivan@119: %
ivan@119: %Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
ivan@119: %quality assessment: From error measurement to structural similarity"
ivan@119: %IEEE Transactios on Image Processing, vol. 13, no. 1, Jan. 2004.
ivan@119: %
ivan@119: %Kindly report any suggestions or corrections to zhouwang@ieee.org
ivan@119: %
ivan@119: %----------------------------------------------------------------------
ivan@119: %
ivan@119: %Input : (1) img1: the first image being compared
ivan@119: %        (2) img2: the second image being compared
ivan@119: %        (3) K: constants in the SSIM index formula (see the above
ivan@119: %            reference). defualt value: K = [0.01 0.03]
ivan@119: %        (4) window: local window for statistics (see the above
ivan@119: %            reference). default widnow is Gaussian given by
ivan@119: %            window = fspecial('gaussian', 11, 1.5);
ivan@119: %        (5) L: dynamic range of the images. default: L = 255
ivan@119: %
ivan@119: %Output: (1) mssim: the mean SSIM index value between 2 images.
ivan@119: %            If one of the images being compared is regarded as 
ivan@119: %            perfect quality, then mssim can be considered as the
ivan@119: %            quality measure of the other image.
ivan@119: %            If img1 = img2, then mssim = 1.
ivan@119: %        (2) ssim_map: the SSIM index map of the test image. The map
ivan@119: %            has a smaller size than the input images. The actual size:
ivan@119: %            size(img1) - size(window) + 1.
ivan@119: %
ivan@119: %Default Usage:
ivan@119: %   Given 2 test images img1 and img2, whose dynamic range is 0-255
ivan@119: %
ivan@119: %   [mssim ssim_map] = ssim_index(img1, img2);
ivan@119: %
ivan@119: %Advanced Usage:
ivan@119: %   User defined parameters. For example
ivan@119: %
ivan@119: %   K = [0.05 0.05];
ivan@119: %   window = ones(8);
ivan@119: %   L = 100;
ivan@119: %   [mssim ssim_map] = ssim_index(img1, img2, K, window, L);
ivan@119: %
ivan@119: %See the results:
ivan@119: %
ivan@119: %   mssim                        %Gives the mssim value
ivan@119: %   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map
ivan@119: %
ivan@119: 
ivan@128: % Removed dependence on Image Processing Toolbox - 
ivan@119: %   Centre for Digital Music, Queen Mary, University of London.
ivan@119: %   This file copyright 2011 Ivan Damnjanovic.
ivan@119: %========================================================================
ivan@119: 
ivan@119: 
ivan@119: if (nargin < 2 || nargin > 5)
ivan@119:    ssim_index = -Inf;
ivan@119:    ssim_map = -Inf;
ivan@119:    return;
ivan@119: end
ivan@119: 
ivan@119: if (size(img1) ~= size(img2))
ivan@119:    ssim_index = -Inf;
ivan@119:    ssim_map = -Inf;
ivan@119:    return;
ivan@119: end
ivan@119: 
ivan@119: [M N] = size(img1);
ivan@119: 
ivan@119: if (nargin == 2)
ivan@119:    if ((M < 11) || (N < 11))
ivan@119: 	   ssim_index = -Inf;
ivan@119: 	   ssim_map = -Inf;
ivan@119:       return
ivan@119:    end
ivan@119: %  window = fspecial('gaussian', 11, 1.5);	%
ivan@119: %
ivan@119: %   Image Processing toolboxdependency 
ivan@119: %   - fspecial function is from IMP toolbox
ivan@119: %   2 - ways to do
ivan@119: % a)
ivan@119: %    [wx,wy] = meshgrid(-5:5, -5:5);
ivan@119: %    window = exp(-(wx.*wx + wy.*wy)/4.5);
ivan@119: %    window = window/sum(window(:));
ivan@119: % b)
ivan@119:    window = [1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06;...
ivan@119:                 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;...
ivan@119:                 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;...
ivan@119:                 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119:                 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119:                 0.000273561160085806,0.00202135875836257,0.00957662749024029,0.0290912256485504,0.0566619704916846,0.0707622377639470,0.0566619704916846,0.0290912256485504,0.00957662749024029,0.00202135875836257,0.000273561160085806;...
ivan@119:                 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119:                 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119:                 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;...
ivan@119:                 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;...
ivan@119:                 1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06];
ivan@119:    K(1) = 0.01;								      % default settings
ivan@119:    K(2) = 0.03;								      %
ivan@119:    L = 255;                                  %
ivan@119: end
ivan@119: 
ivan@119: if (nargin == 3)
ivan@119:    if ((M < 11) || (N < 11))
ivan@119: 	   ssim_index = -Inf;
ivan@119: 	   ssim_map = -Inf;
ivan@119:       return
ivan@119:    end
ivan@119: %  window = fspecial('gaussian', 11, 1.5);	%
ivan@119: %
ivan@119: %   Image Processing toolboxdependency 
ivan@119: %   - fspecial function is from IMP toolbox
ivan@119: %   2 - ways to do
ivan@119: % a)
ivan@119: %    [wx,wy] = meshgrid(-5:5, -5:5);
ivan@119: %    window = exp(-(wx.*wx + wy.*wy)/4.5);
ivan@119: %    window = window/sum(window(:));
ivan@119: % b)
ivan@119:    window = [1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06;...
ivan@119:                 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;...
ivan@119:                 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;...
ivan@119:                 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119:                 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119:                 0.000273561160085806,0.00202135875836257,0.00957662749024029,0.0290912256485504,0.0566619704916846,0.0707622377639470,0.0566619704916846,0.0290912256485504,0.00957662749024029,0.00202135875836257,0.000273561160085806;...
ivan@119:                 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119:                 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119:                 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;...
ivan@119:                 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;...
ivan@119:                 1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06];
ivan@119:    L = 255;
ivan@119:    if (length(K) == 2)
ivan@119:       if (K(1) < 0 || K(2) < 0)
ivan@119: 		   ssim_index = -Inf;
ivan@119:    		ssim_map = -Inf;
ivan@119: 	   	return;
ivan@119:       end
ivan@119:    else
ivan@119: 	   ssim_index = -Inf;
ivan@119:    	ssim_map = -Inf;
ivan@119: 	   return;
ivan@119:    end
ivan@119: end
ivan@119: 
ivan@119: if (nargin == 4)
ivan@119:    [H W] = size(window);
ivan@119:    if ((H*W) < 4 || (H > M) || (W > N))
ivan@119: 	   ssim_index = -Inf;
ivan@119: 	   ssim_map = -Inf;
ivan@119:       return
ivan@119:    end
ivan@119:    L = 255;
ivan@119:    if (length(K) == 2)
ivan@119:       if (K(1) < 0 || K(2) < 0)
ivan@119: 		   ssim_index = -Inf;
ivan@119:    		ssim_map = -Inf;
ivan@119: 	   	return;
ivan@119:       end
ivan@119:    else
ivan@119: 	   ssim_index = -Inf;
ivan@119:    	ssim_map = -Inf;
ivan@119: 	   return;
ivan@119:    end
ivan@119: end
ivan@119: 
ivan@119: if (nargin == 5)
ivan@119:    [H W] = size(window);
ivan@119:    if ((H*W) < 4 || (H > M) || (W > N))
ivan@119: 	   ssim_index = -Inf;
ivan@119: 	   ssim_map = -Inf;
ivan@119:       return
ivan@119:    end
ivan@119:    if (length(K) == 2)
ivan@119:       if (K(1) < 0 || K(2) < 0)
ivan@119: 		   ssim_index = -Inf;
ivan@119:    		ssim_map = -Inf;
ivan@119: 	   	return;
ivan@119:       end
ivan@119:    else
ivan@119: 	   ssim_index = -Inf;
ivan@119:    	ssim_map = -Inf;
ivan@119: 	   return;
ivan@119:    end
ivan@119: end
ivan@119: 
ivan@119: C1 = (K(1)*L)^2;
ivan@119: C2 = (K(2)*L)^2;
ivan@119: window = window/sum(sum(window));
ivan@119: img1 = double(img1);
ivan@119: img2 = double(img2);
ivan@119: 
ivan@119: mu1   = filter2(window, img1, 'valid');
ivan@119: mu2   = filter2(window, img2, 'valid');
ivan@119: mu1_sq = mu1.*mu1;
ivan@119: mu2_sq = mu2.*mu2;
ivan@119: mu1_mu2 = mu1.*mu2;
ivan@119: sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
ivan@119: sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
ivan@119: sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;
ivan@119: 
ivan@119: if (C1 > 0 & C2 > 0)
ivan@119:    ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
ivan@119: else
ivan@119:    numerator1 = 2*mu1_mu2 + C1;
ivan@119:    numerator2 = 2*sigma12 + C2;
ivan@119: 	denominator1 = mu1_sq + mu2_sq + C1;
ivan@119:    denominator2 = sigma1_sq + sigma2_sq + C2;
ivan@119:    ssim_map = ones(size(mu1));
ivan@119:    index = (denominator1.*denominator2 > 0);
ivan@119:    ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
ivan@119:    index = (denominator1 ~= 0) & (denominator2 == 0);
ivan@119:    ssim_map(index) = numerator1(index)./denominator1(index);
ivan@119: end
ivan@119: 
ivan@119: mssim = sum(ssim_map(:), [], 'double')/numel(ssim_map);
ivan@119: 
ivan@119: return