SSIM(structural similarity index) ---图像质量评价指标之结构相似性

Zhou Wang, Alan C. Bovik, Hamid R. Sheikh, and Eero P. Simoncelli. Image Quality Assessment: From Error Visibility to Structural Similarity. IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 13, NO. 4, APRIL 2004

SSIM

SSIM反映了两幅图像之间的相似性，取值为[0,1]。SSIM数值越高越相似。

两幅图像x和y的SSIM计算方法：

其中

K1<<1K_1<<1K1<<1

通常情况下，取α=β=γ=1\alpha=\beta=\gamma=1α=β=γ=1,C3=C2/2C_3= C_2/2C3=C2/2,

此时有

通常K1=0.01,K2=0.03K_1=0.01,K_2=0.03K1=0.01,K2=0.03。

在实际过程中，一般不是对整幅图像求均值、方差、协方差再来计算SSIM。而是对图像分块，滑动窗口 (让stride=1) 求出各窗口的SSIM再来取平均。

除此之外，每个窗口内的均值是利用高斯核函数来求加权均值 ux′u_x^{'}ux′，
方差和协方差都是根据这个加权均值ux′u_x^{'}ux′来求，并且每个元素也需要加权。即 var(x)=∑wi(xi−ux′)2var(x)=\sum w_i(x_i -u_x^{'})^2var(x)=∑wi(xi−ux′)2。

注意每个窗口内生成的卷积模板的权值和都是1。

MATLAB 实现

找到两版认可度较高的matlab代码,其中第二版来自matlab源码:

Code1: 适用于灰度图像

function [mssim, ssim_map] = ssim_index(img1, img2, K, window, L)%========================================================================
%SSIM Index, Version 1.0
%Copyright(c) 2003 Zhou Wang
%All Rights Reserved.
%
%The author is with Howard Hughes Medical Institute, and Laboratory
%for Computational Vision at Center for Neural Science and Courant
%Institute of Mathematical Sciences, New York University.
%
%----------------------------------------------------------------------
%Permission to use, copy, or modify this software and its documentation
%for educational and research purposes only and without fee is hereby
%granted, provided that this copyright notice and the original authors'
%names appear on all copies and supporting documentation. This program
%shall not be used, rewritten, or adapted as the basis of a commercial
%software or hardware product without first obtaining permission of the
%authors. The authors make no representations about the suitability of
%this software for any purpose. It is provided "as is" without express
%or implied warranty.
%----------------------------------------------------------------------
%
%This is an implementation of the algorithm for calculating the
%Structural SIMilarity (SSIM) index between two images. Please refer
%to the following paper:
%
%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
%quality assessment: From error measurement to structural similarity"
%IEEE Transactios on Image Processing, vol. 13, no. 1, Jan. 2004.
%
%Kindly report any suggestions or corrections to zhouwang@ieee.org
%
%----------------------------------------------------------------------
%
%Input : (1) img1: the first image being compared
%        (2) img2: the second image being compared
%        (3) K: constants in the SSIM index formula (see the above
%            reference). defualt value: K = [0.01 0.03]
%        (4) window: local window for statistics (see the above
%            reference). default widnow is Gaussian given by
%            window = fspecial('gaussian', 11, 1.5);
%        (5) L: dynamic range of the images. default: L = 255
%
%Output: (1) mssim: the mean SSIM index value between 2 images.
%            If one of the images being compared is regarded as
%            perfect quality, then mssim can be considered as the
%            quality measure of the other image.
%            If img1 = img2, then mssim = 1.
%        (2) ssim_map: the SSIM index map of the test image. The map
%            has a smaller size than the input images. The actual size:
%            size(img1) - size(window) + 1.
%
%Default Usage:
%   Given 2 test images img1 and img2, whose dynamic range is 0-255
%
%   [mssim ssim_map] = ssim_index(img1, img2);
%
%Advanced Usage:
%   User defined parameters. For example
%
%   K = [0.05 0.05];
%   window = ones(8);                   8\times 8 矩阵
%   L = 100;
%   [mssim ssim_map] = ssim_index(img1, img2, K, window, L);
%
%See the results:
%
%   mssim                        %Gives the mssim value
%   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map
%
%========================================================================if (nargin < 2 || nargin > 5)ssim_index = -Inf;ssim_map = -Inf;return;
endif (size(img1) ~= size(img2))ssim_index = -Inf;ssim_map = -Inf;return;
end[M N] = size(img1);if (nargin == 2)if ((M < 11) || (N < 11))ssim_index = -Inf;ssim_map = -Inf;returnendwindow = fspecial('gaussian', 11, 1.5); %K(1) = 0.01;                                    % default settingsK(2) = 0.03;                                     %L = 255;                                  %
endif (nargin == 3)if ((M < 11) || (N < 11))ssim_index = -Inf;ssim_map = -Inf;returnendwindow = fspecial('gaussian', 11, 1.5);  % 11 \times 11 的矩阵L = 255;if (length(K) == 2)if (K(1) < 0 || K(2) < 0)ssim_index = -Inf;ssim_map = -Inf;return;endelsessim_index = -Inf;ssim_map = -Inf;return;end
endif (nargin == 4)[H W] = size(window);if ((H*W) < 4 || (H > M) || (W > N))ssim_index = -Inf;ssim_map = -Inf;returnendL = 255;if (length(K) == 2)if (K(1) < 0 || K(2) < 0)ssim_index = -Inf;ssim_map = -Inf;return;endelsessim_index = -Inf;ssim_map = -Inf;return;end
endif (nargin == 5)[H W] = size(window);if ((H*W) < 4 || (H > M) || (W > N))ssim_index = -Inf;ssim_map = -Inf;returnendif (length(K) == 2)if (K(1) < 0 || K(2) < 0)ssim_index = -Inf;ssim_map = -Inf;return;endelsessim_index = -Inf;ssim_map = -Inf;return;end
endC1 = (K(1)*L)^2;
C2 = (K(2)*L)^2;
window = window/sum(sum(window));                      %使权值和为1
img1 = double(img1);
img2 = double(img2);mu1   = filter2(window, img1, 'valid');     % [M-filter_size+1 N-filter_size+1 ] 大小
mu2   = filter2(window, img2, 'valid');     % [M-filter_size+1 N-filter_size+1 ] 大小
mu1_sq = mu1.*mu1;
mu2_sq = mu2.*mu2;
mu1_mu2 = mu1.*mu2;
sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;
sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;
sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;if (C1 > 0 & C2 > 0)ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
elsenumerator1 = 2*mu1_mu2 + C1;numerator2 = 2*sigma12 + C2;denominator1 = mu1_sq + mu2_sq + C1;denominator2 = sigma1_sq + sigma2_sq + C2;ssim_map = ones(size(mu1));index = (denominator1.*denominator2 > 0);ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));index = (denominator1 ~= 0) & (denominator2 == 0);ssim_map(index) = numerator1(index)./denominator1(index);
endmssim = mean2(ssim_map);return

第二份代码之所以可以处理3通道图像的原因是生成了一个三维的高斯核函数卷积模板。

Code2: 适用于2通道或3通道图像

function [ssimval, ssimmap] = ssim(varargin)
%SSIM Structural Similarity Index for measuring image quality
%   SSIMVAL = SSIM(A, REF) calculates the Structural Similarity Index
%   (SSIM) value for image A, with the image REF as the reference. A and
%   REF can be 2D grayscale or 3D volume images, and must be of the same
%   size and class.
%
%   [SSIMVAL, SSIMMAP] = SSIM(A, REF) also returns the local SSIM value for
%   each pixel in SSIMMAP. SSIMMAP has the same size as A.
%
%   [SSIMVAL, SSIMMAP] = SSIM(A, REF, NAME1, VAL1,...) calculates the SSIM
%   value using name-value pairs to control aspects of the computation.
%   Parameter names can be abbreviated.
%
%   Parameters include:
%
%   'Radius'                 - Specifies the standard deviation of
%                              isotropic Gaussian function used for
%                              weighting the neighborhood pixels around a
%                              pixel for estimating local statistics. This
%                              weighting is used to avoid blocking
%                              artifacts in estimating local statistics.
%                              The default value is 1.5.
%
%   'DynamicRange'           - Positive scalar, L, that specifies the
%                              dynamic range of the input image. By
%                              default, L is chosen based on the class of
%                              the input image A, as L =
%                              diff(getrangefromclass(A)). Note that when
%                              class of A is single or double, L = 1 by
%                              default.
%
%   'RegularizationConstants'- Three-element vector, [C1 C2 C3], of
%                              non-negative real numbers that specifies the
%                              regularization constants for the luminance,
%                              contrast, and structural terms (see [1]),
%                              respectively. The regularization constants
%                              are used to avoid instability for image
%                              regions where the local mean or standard
%                              deviation is close to zero. Therefore, small
%                              non-zero values should be used for these
%                              constants. By default, C1 = (0.01*L).^2, C2
%                              = (0.03*L).^2, and C3 = C2/2, where L is the
%                              specified 'DynamicRange' value. If a value
%                              of 'DynamicRange' is not specified, the
%                              default value is used (see name-value pair
%                              'DynamicRange').
%
%   'Exponents'               - Three-element vector [alpha beta gamma],
%                               of non-negative real numbers that specifies
%                               the exponents for the luminance, contrast,
%                               and structural terms (see [1]),
%                               respectively. By default, all the three
%                               exponents are 1, i.e. the vector is [1 1
%                               1].
%
%   Notes
%   -----
%   1. A and REF can be arrays of upto three dimensions. All 3D arrays
%      are considered 3D volumetric images. RGB images will also be
%      processed as 3D volumetric images.
%
%   2. Input image A and reference image REF are converted to
%      floating-point type for internal computation.
%
%   3. For signed-integer images (int16), an offset is applied to bring the
%      gray values in the non-negative range before computing the SSIM
%      index.
%
%   Example
%   ---------
%   This example shows how to compute SSIM value for a blurred image given
%   the original reference image.
%
%   ref = imread('pout.tif');
%   H = fspecial('Gaussian',[11 11],1.5);
%   A = imfilter(ref,H,'replicate');  % 'replicate' 不改变size
%
%   subplot(1,2,1); imshow(ref); title('Reference Image');
%   subplot(1,2,2); imshow(A);   title('Blurred Image');
%
%   [ssimval, ssimmap] = ssim(A,ref);
%
%   fprintf('The SSIM value is %0.4f.\n',ssimval);
%
%   figure, imshow(ssimmap,[]);
%   title(sprintf('SSIM Index Map - Mean SSIM Value is %0.4f',ssimval));
%
%   Class Support
%   -------------
%   Input arrays A and REF must be one of the following classes: uint8,
%   int16, uint16, single, or double. Both A and REF must be of the same
%   class. They must be nonsparse. SSIMVAL is a scalar and SSIMMAP is an
%   array of the same size as A. Both SSIMVAL and SSIMMAP are of class
%   double, unless A and REF are of class single in which case SSIMVAL and
%   SSIMMAP are of class single.
%
%   References:
%   -----------
%   [1] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
%       Quality Assessment: From Error Visibility to Structural
%       Similarity," IEEE Transactions on Image Processing, Volume 13,
%       Issue 4, pp. 600- 612, 2004.
%
%   See also IMMSE, MEAN, MEDIAN, PSNR, SUM, VAR.%   Copyright 2013-2017 The MathWorks, Inc.narginchk(2,10);
%narginchk(minArgs,maxArgs) 验证当前执行的函数调用中的输入参数数目。如果调用中指定的输入数目小于 minArgs 或大于 maxArgs，narginchk 将引发错误。
%如果输入数目在 minArgs 与 maxArgs 之间（包括二者），则 narginchk 不会执行任何操作。args = matlab.images.internal.stringToChar(varargin);
[A, ref, C, exponents, radius] = parse_inputs(args{:});if isempty(A)ssimval = zeros(0, 'like', A);ssimmap = A;return;
endif isa(A,'int16') % int16 is the only allowed signed-integer type for A and ref.% Add offset for signed-integer types to bring values in the% non-negative range.A = double(A) - double(intmin('int16'));ref = double(ref) - double(intmin('int16'));
elseif isinteger(A)A = double(A);ref = double(ref);
end% Gaussian weighting function
gaussFilt = getGaussianWeightingFilter(radius,ndims(A));% Weighted-mean and weighted-variance computations
mux2 = imfilter(A, gaussFilt,'conv','replicate');
muy2 = imfilter(ref, gaussFilt,'conv','replicate');
muxy = mux2.*muy2;
mux2 = mux2.^2;
muy2 = muy2.^2;sigmax2 = imfilter(A.^2,gaussFilt,'conv','replicate') - mux2;
sigmay2 = imfilter(ref.^2,gaussFilt,'conv','replicate') - muy2;
sigmaxy = imfilter(A.*ref,gaussFilt,'conv','replicate') - muxy;% Compute SSIM index
if (C(3) == C(2)/2) && isequal(exponents(:),ones(3,1))% Special case: Equation 13 from [1]num = (2*muxy + C(1)).*(2*sigmaxy + C(2));den = (mux2 + muy2 + C(1)).*(sigmax2 + sigmay2 + C(2));if (C(1) > 0) && (C(2) > 0)ssimmap = num./den;else% Need to guard against divide-by-zero if either C(1) or C(2) is 0.isDenNonZero = (den ~= 0);           ssimmap = ones(size(A));ssimmap(isDenNonZero) = num(isDenNonZero)./den(isDenNonZero);endelse% General case: Equation 12 from [1] % Luminance termif (exponents(1) > 0)num = 2*muxy + C(1);den = mux2 + muy2 + C(1); ssimmap = guardedDivideAndExponent(num,den,C(1),exponents(1));else ssimmap = ones(size(A), 'like', A);end% Contrast termsigmaxsigmay = [];if (exponents(2) > 0)  sigmaxsigmay = sqrt(sigmax2.*sigmay2);num = 2*sigmaxsigmay + C(2);den = sigmax2 + sigmay2 + C(2); ssimmap = ssimmap.*guardedDivideAndExponent(num,den,C(2),exponents(2));        end% Structure termif (exponents(3) > 0)num = sigmaxy + C(3);if isempty(sigmaxsigmay)sigmaxsigmay = sqrt(sigmax2.*sigmay2);endden = sigmaxsigmay + C(3); ssimmap = ssimmap.*guardedDivideAndExponent(num,den,C(3),exponents(3));        endendssimval = mean(ssimmap(:));end% -------------------------------------------------------------------------
function component = guardedDivideAndExponent(num, den, C, exponent)if C > 0component = num./den;
elsecomponent = ones(size(num),'like',num);isDenNonZero = (den ~= 0);component(isDenNonZero) = num(isDenNonZero)./den(isDenNonZero);
endif (exponent ~= 1)component = component.^exponent;
endendfunction gaussFilt = getGaussianWeightingFilter(radius,N)
% Get 2D or 3D Gaussian weighting filterfiltRadius = ceil(radius*3); % 3 Standard deviations include >99% of the area.
filtSize = 2*filtRadius + 1;if (N < 3)% 2D Gaussian mask can be used for filtering even one-dimensional% signals using imfilter. gaussFilt = fspecial('gaussian',[filtSize filtSize],radius);
else % 3D Gaussian mask[x,y,z] = ndgrid(-filtRadius:filtRadius,-filtRadius:filtRadius, ...-filtRadius:filtRadius);arg = -(x.*x + y.*y + z.*z)/(2*radius*radius);gaussFilt = exp(arg);gaussFilt(gaussFilt<eps*max(gaussFilt(:))) = 0;sumFilt = sum(gaussFilt(:));if (sumFilt ~= 0)gaussFilt  = gaussFilt/sumFilt;end
%总之返回的系数模板权值和为1
endendfunction [A, ref, C, exponents, radius] = parse_inputs(varargin)validImageTypes = {'uint8','uint16','int16','single','double'};A = varargin{1};
validateattributes(A,validImageTypes,{'nonsparse','real'},mfilename,'A',1);ref = varargin{2};
validateattributes(ref,validImageTypes,{'nonsparse','real'},mfilename,'REF',2);if ~isa(A,class(ref))error(message('images:validate:differentClassMatrices','A','REF'));
endif ~isequal(size(A),size(ref))error(message('images:validate:unequalSizeMatrices','A','REF'));
endif (ndims(A) > 3)error(message('images:validate:tooManyDimensions','A and REF',3));
end% Default values for parameters
dynmRange = diff(getrangefromclass(A));
C = [];
exponents = [1 1 1];
radius = 1.5;args_names = {'dynamicrange', 'regularizationconstants','exponents',...'radius'};for i = 3:2:narginarg = varargin{i};if ischar(arg)        idx = find(strncmpi(arg, args_names, numel(arg)));if isempty(idx)error(message('images:validate:unknownInputString', arg))elseif numel(idx) > 1error(message('images:validate:ambiguousInputString', arg))elseif numel(idx) == 1if (i+1 > nargin) error(message('images:validate:missingParameterValue'));             endif idx == 1dynmRange = varargin{i+1};validateattributes(dynmRange,{'numeric'},{'positive', ...'finite', 'real', 'nonempty','scalar'}, mfilename, ...'DynamicRange',i);dynmRange = double(dynmRange);elseif idx == 2C = varargin{i+1};validateattributes(C,{'numeric'},{'nonnegative','finite', ...'real','nonempty','vector', 'numel', 3}, mfilename, ...'RegularizationConstants',i);                              C = double(C);                              elseif idx == 3exponents = varargin{i+1};validateattributes(exponents,{'numeric'},{'nonnegative', ...'finite', 'real', 'nonempty','vector', 'numel', 3}, ...mfilename,'Exponents',i);exponents = double(exponents);elseif idx == 4radius = varargin{i+1};validateattributes(radius,{'numeric'},{'positive','finite', ...'real', 'nonempty','scalar'}, mfilename,'Radius',i);radius = double(radius);endend    elseerror(message('images:validate:mustBeString')); end
end% If 'RegularizationConstants' is not specified, choose default C.
if isempty(C)C = [(0.01*dynmRange).^2 (0.03*dynmRange).^2 ((0.03*dynmRange).^2)/2];
endend

对于代码核心部分的理解

mux2 = imfilter(A, gaussFilt,'conv','replicate');
muy2 = imfilter(ref, gaussFilt,'conv','replicate');
muxy = mux2.*muy2;
mux2 = mux2.^2;
muy2 = muy2.^2;sigmax2 = imfilter(A.^2,gaussFilt,'conv','replicate') - mux2;
sigmay2 = imfilter(ref.^2,gaussFilt,'conv','replicate') - muy2;
sigmaxy = imfilter(A.*ref,gaussFilt,'conv','replicate') - muxy;

ux′u_x'ux′为加权均值,

var(x)=∑wi(xi−ux′)2=∑wi(xi2−2ux′xi+ux′2)=wixi2−2ux′∑wixi+∑wiux′2=wixi2+ux′(ux′∑wi−2∑wixi)=∑wixi2+ux′(ux′−2ux′)=∑wixi2−ux′2\begin{aligned} var(x) &= \sum w_i(x_i-u_x^{'})^2 \\ &= \sum w_i (x_i^2-2u_x^{'}x_i +u_x^{'2}) \\ &= w_i x_i^2 -2u_x^{'} \sum w_i x_i +\sum w_i u_x^{'2} \\ &= w_i x_i^2 + u_x^{'} (u_x^{'} \sum w_i- 2\sum w_i x_i )\\ &= \sum w_i x_i^2 + u_x^{'} (u_x^{'} -2 u_x^{'}) \\ &= \sum w_i x_i^2 - u_x^{'2} \end{aligned} var(x)=∑wi(xi−ux′)2=∑wi(xi2−2ux′xi+ux′2)=wixi2−2ux′∑wixi+∑wiux′2=wixi2+ux′(ux′∑wi−2∑wixi)=∑wixi2+ux′(ux′−2ux′)=∑wixi2−ux′2

注意∑wi=1\sum w_i=1∑wi=1。

cov(x,y)=∑wi(xi−ux′)(yi−uy′)=∑wi(xiyi−xiuy′−ux′yi+ux′uy′)=∑wixiyi+ux′uy′−uy′∑wixi−ux′∑wiyi=∑wixiyi+ux′uy′−uy′ux′−ux′uy′=∑wixiyi−ux′uy′\begin{aligned} cov(x,y) &= \sum w_i(x_i-u_x^{'}) (y_i-u_y^{'}) \\ &= \sum w_i (x_i y_i - x_i u_y^{'} - u_x^{'} y_i+ u_x^{'} u_y^{'}) \\ &= \sum w_i x_i y_i + u_x^{'} u_y^{'}-u_y^{'} \sum w_i x_i-u_x^{'} \sum w_i y_i \\ &= \sum w_i x_i y_i + u_x^{'} u_y^{'}-u_y^{'} u_x^{'} - u_x^{'} u_y^{'}\\ &= \sum w_i x_i y_i - u_x^{'} u_y^{'}\\ \end{aligned} cov(x,y)=∑wi(xi−ux′)(yi−uy′)=∑wi(xiyi−xiuy′−ux′yi+ux′uy′)=∑wixiyi+ux′uy′−uy′∑wixi−ux′∑wiyi=∑wixiyi+ux′uy′−uy′ux′−ux′uy′=∑wixiyi−ux′uy′