THESE/codes/ssim.m

   1 function [mssim, ssim_map] = ssim(img1, img2, K, window, L)\r
   2 \r
   3 % ========================================================================\r
   4 % SSIM Index with automatic downsampling, Version 1.0\r
   5 % Copyright(c) 2009 Zhou Wang\r
   6 % All Rights Reserved.\r
   7 %\r
   8 % ----------------------------------------------------------------------\r
   9 % Permission to use, copy, or modify this software and its documentation\r
  10 % for educational and research purposes only and without fee is hereby\r
  11 % granted, provided that this copyright notice and the original authors'\r
  12 % names appear on all copies and supporting documentation. This program\r
  13 % shall not be used, rewritten, or adapted as the basis of a commercial\r
  14 % software or hardware product without first obtaining permission of the\r
  15 % authors. The authors make no representations about the suitability of\r
  16 % this software for any purpose. It is provided "as is" without express\r
  17 % or implied warranty.\r
  18 %----------------------------------------------------------------------\r
  19 %\r
  20 % This is an implementation of the algorithm for calculating the\r
  21 % Structural SIMilarity (SSIM) index between two images\r
  22 %\r
  23 % Please refer to the following paper and the website with suggested usage\r
  24 %\r
  25 % Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image\r
  26 % quality assessment: From error visibility to structural similarity,"\r
  27 % IEEE Transactios on Image Processing, vol. 13, no. 4, pp. 600-612,\r
  28 % Apr. 2004.\r
  29 %\r
  30 % http://www.ece.uwaterloo.ca/~z70wang/research/ssim/\r
  31 %\r
  32 % Note: This program is different from ssim_index.m, where no automatic\r
  33 % downsampling is performed. (downsampling was done in the above paper\r
  34 % and was described as suggested usage in the above website.)\r
  35 %\r
  36 % Kindly report any suggestions or corrections to zhouwang@ieee.org\r
  37 %\r
  38 %----------------------------------------------------------------------\r
  39 %\r
  40 %Input : (1) img1: the first image being compared\r
  41 %        (2) img2: the second image being compared\r
  42 %        (3) K: constants in the SSIM index formula (see the above\r
  43 %            reference). defualt value: K = [0.01 0.03]\r
  44 %        (4) window: local window for statistics (see the above\r
  45 %            reference). default widnow is Gaussian given by\r
  46 %            window = fspecial('gaussian', 11, 1.5);\r
  47 %        (5) L: dynamic range of the images. default: L = 255\r
  48 %\r
  49 %Output: (1) mssim: the mean SSIM index value between 2 images.\r
  50 %            If one of the images being compared is regarded as \r
  51 %            perfect quality, then mssim can be considered as the\r
  52 %            quality measure of the other image.\r
  53 %            If img1 = img2, then mssim = 1.\r
  54 %        (2) ssim_map: the SSIM index map of the test image. The map\r
  55 %            has a smaller size than the input images. The actual size\r
  56 %            depends on the window size and the downsampling factor.\r
  57 %\r
  58 %Basic Usage:\r
  59 %   Given 2 test images img1 and img2, whose dynamic range is 0-255\r
  60 %\r
  61 %   [mssim, ssim_map] = ssim(img1, img2);\r
  62 %\r
  63 %Advanced Usage:\r
  64 %   User defined parameters. For example\r
  65 %\r
  66 %   K = [0.05 0.05];\r
  67 %   window = ones(8);\r
  68 %   L = 100;\r
  69 %   [mssim, ssim_map] = ssim(img1, img2, K, window, L);\r
  70 %\r
  71 %Visualize the results:\r
  72 %\r
  73 %   mssim                        %Gives the mssim value\r
  74 %   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map\r
  75 %========================================================================\r
  76 \r
  77 \r
  78 if (nargin < 2 || nargin > 5)\r
  79    mssim = -Inf;\r
  80    ssim_map = -Inf;\r
  81    return;\r
  82 end\r
  83 \r
  84 if (size(img1) ~= size(img2))\r
  85    mssim = -Inf;\r
  86    ssim_map = -Inf;\r
  87    return;\r
  88 end\r
  89 \r
  90 [M N] = size(img1);\r
  91 \r
  92 if (nargin == 2)\r
  93    if ((M < 11) || (N < 11))\r
  94            mssim = -Inf;\r
  95            ssim_map = -Inf;\r
  96       return\r
  97    end\r
  98    window = fspecial('gaussian', 11, 1.5);      %\r
  99    K(1) = 0.01;                                 % default settings\r
 100    K(2) = 0.03;                                 %\r
 101    L = 255;                                     %\r
 102 end\r
 103 \r
 104 if (nargin == 3)\r
 105    if ((M < 11) || (N < 11))\r
 106            mssim = -Inf;\r
 107            ssim_map = -Inf;\r
 108       return\r
 109    end\r
 110    window = fspecial('gaussian', 11, 1.5);\r
 111    L = 255;\r
 112    if (length(K) == 2)\r
 113       if (K(1) < 0 || K(2) < 0)\r
 114                    mssim = -Inf;\r
 115                 ssim_map = -Inf;\r
 116                 return;\r
 117       end\r
 118    else\r
 119            mssim = -Inf;\r
 120         ssim_map = -Inf;\r
 121            return;\r
 122    end\r
 123 end\r
 124 \r
 125 if (nargin == 4)\r
 126    [H W] = size(window);\r
 127    if ((H*W) < 4 || (H > M) || (W > N))\r
 128            mssim = -Inf;\r
 129            ssim_map = -Inf;\r
 130       return\r
 131    end\r
 132    L = 255;\r
 133    if (length(K) == 2)\r
 134       if (K(1) < 0 || K(2) < 0)\r
 135                    mssim = -Inf;\r
 136                 ssim_map = -Inf;\r
 137                 return;\r
 138       end\r
 139    else\r
 140            mssim = -Inf;\r
 141         ssim_map = -Inf;\r
 142            return;\r
 143    end\r
 144 end\r
 145 \r
 146 if (nargin == 5)\r
 147    [H W] = size(window);\r
 148    if ((H*W) < 4 || (H > M) || (W > N))\r
 149            mssim = -Inf;\r
 150            ssim_map = -Inf;\r
 151       return\r
 152    end\r
 153    if (length(K) == 2)\r
 154       if (K(1) < 0 || K(2) < 0)\r
 155                    mssim = -Inf;\r
 156                 ssim_map = -Inf;\r
 157                 return;\r
 158       end\r
 159    else\r
 160            mssim = -Inf;\r
 161         ssim_map = -Inf;\r
 162            return;\r
 163    end\r
 164 end\r
 165 \r
 166 \r
 167 img1 = double(img1);\r
 168 img2 = double(img2);\r
 169 \r
 170 % automatic downsampling\r
 171 f = max(1,round(min(M,N)/256));\r
 172 %downsampling by f\r
 173 %use a simple low-pass filter \r
 174 if(f>1)\r
 175     lpf = ones(f,f);\r
 176     lpf = lpf/sum(lpf(:));\r
 177     img1 = imfilter(img1,lpf,'symmetric','same');\r
 178     img2 = imfilter(img2,lpf,'symmetric','same');\r
 179 \r
 180     img1 = img1(1:f:end,1:f:end);\r
 181     img2 = img2(1:f:end,1:f:end);\r
 182 end\r
 183 \r
 184 C1 = (K(1)*L)^2;\r
 185 C2 = (K(2)*L)^2;\r
 186 window = window/sum(sum(window));\r
 187 \r
 188 mu1   = filter2(window, img1, 'valid');\r
 189 mu2   = filter2(window, img2, 'valid');\r
 190 mu1_sq = mu1.*mu1;\r
 191 mu2_sq = mu2.*mu2;\r
 192 mu1_mu2 = mu1.*mu2;\r
 193 sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;\r
 194 sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;\r
 195 sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;\r
 196 \r
 197 if (C1 > 0 && C2 > 0)\r
 198    ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));\r
 199 else\r
 200    numerator1 = 2*mu1_mu2 + C1;\r
 201    numerator2 = 2*sigma12 + C2;\r
 202         denominator1 = mu1_sq + mu2_sq + C1;\r
 203    denominator2 = sigma1_sq + sigma2_sq + C2;\r
 204    ssim_map = ones(size(mu1));\r
 205    index = (denominator1.*denominator2 > 0);\r
 206    ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));\r
 207    index = (denominator1 ~= 0) & (denominator2 == 0);\r
 208    ssim_map(index) = numerator1(index)./denominator1(index);\r
 209 end\r
 210 \r
 211 mssim = mean2(ssim_map);\r
 212 \r
 213 return