comparison util/SMALL_ssim_index.m @ 119:5356a6e13a25 sup_163_ssim_IMG_T_dependeny

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