annotate util/SMALL_ssim_index.m @ 183:0d7a81655ef2 danieleb

removed cumulative coherence calculation
author Daniele Barchiesi <daniele.barchiesi@eecs.qmul.ac.uk>
date Fri, 27 Jan 2012 13:15:11 +0000
parents 8e660fd14774
children
rev   line source
ivan@119 1 function [mssim, ssim_map] = SMALL_ssim_index(img1, img2, K, window, L)
ivan@128 2 %% SSIM Index - Removed dependence on Image Processing Toolbox
ivan@119 3 %========================================================================
ivan@119 4 %SSIM Index, Version 1.0
ivan@119 5 %Copyright(c) 2003 Zhou Wang
ivan@119 6 %All Rights Reserved.
ivan@119 7 %
ivan@119 8 %The author is with Howard Hughes Medical Institute, and Laboratory
ivan@119 9 %for Computational Vision at Center for Neural Science and Courant
ivan@119 10 %Institute of Mathematical Sciences, New York University.
ivan@119 11 %
ivan@119 12 %----------------------------------------------------------------------
ivan@119 13 %Permission to use, copy, or modify this software and its documentation
ivan@119 14 %for educational and research purposes only and without fee is hereby
ivan@119 15 %granted, provided that this copyright notice and the original authors'
ivan@119 16 %names appear on all copies and supporting documentation. This program
ivan@119 17 %shall not be used, rewritten, or adapted as the basis of a commercial
ivan@119 18 %software or hardware product without first obtaining permission of the
ivan@119 19 %authors. The authors make no representations about the suitability of
ivan@119 20 %this software for any purpose. It is provided "as is" without express
ivan@119 21 %or implied warranty.
ivan@119 22 %----------------------------------------------------------------------
ivan@119 23 %
ivan@119 24 %This is an implementation of the algorithm for calculating the
ivan@119 25 %Structural SIMilarity (SSIM) index between two images. Please refer
ivan@119 26 %to the following paper:
ivan@119 27 %
ivan@119 28 %Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
ivan@119 29 %quality assessment: From error measurement to structural similarity"
ivan@119 30 %IEEE Transactios on Image Processing, vol. 13, no. 1, Jan. 2004.
ivan@119 31 %
ivan@119 32 %Kindly report any suggestions or corrections to zhouwang@ieee.org
ivan@119 33 %
ivan@119 34 %----------------------------------------------------------------------
ivan@119 35 %
ivan@119 36 %Input : (1) img1: the first image being compared
ivan@119 37 % (2) img2: the second image being compared
ivan@119 38 % (3) K: constants in the SSIM index formula (see the above
ivan@119 39 % reference). defualt value: K = [0.01 0.03]
ivan@119 40 % (4) window: local window for statistics (see the above
ivan@119 41 % reference). default widnow is Gaussian given by
ivan@119 42 % window = fspecial('gaussian', 11, 1.5);
ivan@119 43 % (5) L: dynamic range of the images. default: L = 255
ivan@119 44 %
ivan@119 45 %Output: (1) mssim: the mean SSIM index value between 2 images.
ivan@119 46 % If one of the images being compared is regarded as
ivan@119 47 % perfect quality, then mssim can be considered as the
ivan@119 48 % quality measure of the other image.
ivan@119 49 % If img1 = img2, then mssim = 1.
ivan@119 50 % (2) ssim_map: the SSIM index map of the test image. The map
ivan@119 51 % has a smaller size than the input images. The actual size:
ivan@119 52 % size(img1) - size(window) + 1.
ivan@119 53 %
ivan@119 54 %Default Usage:
ivan@119 55 % Given 2 test images img1 and img2, whose dynamic range is 0-255
ivan@119 56 %
ivan@119 57 % [mssim ssim_map] = ssim_index(img1, img2);
ivan@119 58 %
ivan@119 59 %Advanced Usage:
ivan@119 60 % User defined parameters. For example
ivan@119 61 %
ivan@119 62 % K = [0.05 0.05];
ivan@119 63 % window = ones(8);
ivan@119 64 % L = 100;
ivan@119 65 % [mssim ssim_map] = ssim_index(img1, img2, K, window, L);
ivan@119 66 %
ivan@119 67 %See the results:
ivan@119 68 %
ivan@119 69 % mssim %Gives the mssim value
ivan@119 70 % imshow(max(0, ssim_map).^4) %Shows the SSIM index map
ivan@119 71 %
ivan@119 72
ivan@128 73 % Removed dependence on Image Processing Toolbox -
ivan@119 74 % Centre for Digital Music, Queen Mary, University of London.
ivan@119 75 % This file copyright 2011 Ivan Damnjanovic.
ivan@119 76 %========================================================================
ivan@119 77
ivan@119 78
ivan@119 79 if (nargin < 2 || nargin > 5)
ivan@119 80 ssim_index = -Inf;
ivan@119 81 ssim_map = -Inf;
ivan@119 82 return;
ivan@119 83 end
ivan@119 84
ivan@119 85 if (size(img1) ~= size(img2))
ivan@119 86 ssim_index = -Inf;
ivan@119 87 ssim_map = -Inf;
ivan@119 88 return;
ivan@119 89 end
ivan@119 90
ivan@119 91 [M N] = size(img1);
ivan@119 92
ivan@119 93 if (nargin == 2)
ivan@119 94 if ((M < 11) || (N < 11))
ivan@119 95 ssim_index = -Inf;
ivan@119 96 ssim_map = -Inf;
ivan@119 97 return
ivan@119 98 end
ivan@119 99 % window = fspecial('gaussian', 11, 1.5); %
ivan@119 100 %
ivan@119 101 % Image Processing toolboxdependency
ivan@119 102 % - fspecial function is from IMP toolbox
ivan@119 103 % 2 - ways to do
ivan@119 104 % a)
ivan@119 105 % [wx,wy] = meshgrid(-5:5, -5:5);
ivan@119 106 % window = exp(-(wx.*wx + wy.*wy)/4.5);
ivan@119 107 % window = window/sum(window(:));
ivan@119 108 % b)
ivan@119 109 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 110 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 111 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 112 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119 113 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119 114 0.000273561160085806,0.00202135875836257,0.00957662749024029,0.0290912256485504,0.0566619704916846,0.0707622377639470,0.0566619704916846,0.0290912256485504,0.00957662749024029,0.00202135875836257,0.000273561160085806;...
ivan@119 115 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119 116 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119 117 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 118 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 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 120 K(1) = 0.01; % default settings
ivan@119 121 K(2) = 0.03; %
ivan@119 122 L = 255; %
ivan@119 123 end
ivan@119 124
ivan@119 125 if (nargin == 3)
ivan@119 126 if ((M < 11) || (N < 11))
ivan@119 127 ssim_index = -Inf;
ivan@119 128 ssim_map = -Inf;
ivan@119 129 return
ivan@119 130 end
ivan@119 131 % window = fspecial('gaussian', 11, 1.5); %
ivan@119 132 %
ivan@119 133 % Image Processing toolboxdependency
ivan@119 134 % - fspecial function is from IMP toolbox
ivan@119 135 % 2 - ways to do
ivan@119 136 % a)
ivan@119 137 % [wx,wy] = meshgrid(-5:5, -5:5);
ivan@119 138 % window = exp(-(wx.*wx + wy.*wy)/4.5);
ivan@119 139 % window = window/sum(window(:));
ivan@119 140 % b)
ivan@119 141 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 142 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 143 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 144 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119 145 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119 146 0.000273561160085806,0.00202135875836257,0.00957662749024029,0.0290912256485504,0.0566619704916846,0.0707622377639470,0.0566619704916846,0.0290912256485504,0.00957662749024029,0.00202135875836257,0.000273561160085806;...
ivan@119 147 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;...
ivan@119 148 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;...
ivan@119 149 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 150 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 151 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 152 L = 255;
ivan@119 153 if (length(K) == 2)
ivan@119 154 if (K(1) < 0 || K(2) < 0)
ivan@119 155 ssim_index = -Inf;
ivan@119 156 ssim_map = -Inf;
ivan@119 157 return;
ivan@119 158 end
ivan@119 159 else
ivan@119 160 ssim_index = -Inf;
ivan@119 161 ssim_map = -Inf;
ivan@119 162 return;
ivan@119 163 end
ivan@119 164 end
ivan@119 165
ivan@119 166 if (nargin == 4)
ivan@119 167 [H W] = size(window);
ivan@119 168 if ((H*W) < 4 || (H > M) || (W > N))
ivan@119 169 ssim_index = -Inf;
ivan@119 170 ssim_map = -Inf;
ivan@119 171 return
ivan@119 172 end
ivan@119 173 L = 255;
ivan@119 174 if (length(K) == 2)
ivan@119 175 if (K(1) < 0 || K(2) < 0)
ivan@119 176 ssim_index = -Inf;
ivan@119 177 ssim_map = -Inf;
ivan@119 178 return;
ivan@119 179 end
ivan@119 180 else
ivan@119 181 ssim_index = -Inf;
ivan@119 182 ssim_map = -Inf;
ivan@119 183 return;
ivan@119 184 end
ivan@119 185 end
ivan@119 186
ivan@119 187 if (nargin == 5)
ivan@119 188 [H W] = size(window);
ivan@119 189 if ((H*W) < 4 || (H > M) || (W > N))
ivan@119 190 ssim_index = -Inf;
ivan@119 191 ssim_map = -Inf;
ivan@119 192 return
ivan@119 193 end
ivan@119 194 if (length(K) == 2)
ivan@119 195 if (K(1) < 0 || K(2) < 0)
ivan@119 196 ssim_index = -Inf;
ivan@119 197 ssim_map = -Inf;
ivan@119 198 return;
ivan@119 199 end
ivan@119 200 else
ivan@119 201 ssim_index = -Inf;
ivan@119 202 ssim_map = -Inf;
ivan@119 203 return;
ivan@119 204 end
ivan@119 205 end
ivan@119 206
ivan@119 207 C1 = (K(1)*L)^2;
ivan@119 208 C2 = (K(2)*L)^2;
ivan@119 209 window = window/sum(sum(window));
ivan@119 210 img1 = double(img1);
ivan@119 211 img2 = double(img2);
ivan@119 212
ivan@119 213 mu1 = filter2(window, img1, 'valid');
ivan@119 214 mu2 = filter2(window, img2, 'valid');
ivan@119 215 mu1_sq = mu1.*mu1;
ivan@119 216 mu2_sq = mu2.*mu2;
ivan@119 217 mu1_mu2 = mu1.*mu2;
ivan@119 218 sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
ivan@119 219 sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
ivan@119 220 sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;
ivan@119 221
ivan@119 222 if (C1 > 0 & C2 > 0)
ivan@119 223 ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
ivan@119 224 else
ivan@119 225 numerator1 = 2*mu1_mu2 + C1;
ivan@119 226 numerator2 = 2*sigma12 + C2;
ivan@119 227 denominator1 = mu1_sq + mu2_sq + C1;
ivan@119 228 denominator2 = sigma1_sq + sigma2_sq + C2;
ivan@119 229 ssim_map = ones(size(mu1));
ivan@119 230 index = (denominator1.*denominator2 > 0);
ivan@119 231 ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
ivan@119 232 index = (denominator1 ~= 0) & (denominator2 == 0);
ivan@119 233 ssim_map(index) = numerator1(index)./denominator1(index);
ivan@119 234 end
ivan@119 235
ivan@119 236 mssim = sum(ssim_map(:), [], 'double')/numel(ssim_map);
ivan@119 237
ivan@119 238 return