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计算方法：
其中
K
1
<
<
1
K_1<<1
K1<<1
通常情况下，取 α = β = γ = 1 \alpha=\beta=\gamma=1 α=β=γ=1, C 3 = C 2 / 2 C_3= C_2/2 C3=C2/2,
此时有
通常 K 1 = 0.01 , K 2 = 0.03 K_1=0.01,K_2=0.03 K1=0.01,K2=0.03。
在实际过程中，一般不是对整幅图像求均值、方差、协方差再来计算SSIM。而是对图像分块，滑动窗口 (让stride=1) 求出各窗口的SSIM再来取平均。
除此之外，每个窗口内的均值是利用高斯核函数来求加权均值
u
x
′
u_x^{'}
ux′，
方差和协方差都是根据这个加权均值
u
x
′
u_x^{'}
ux′来求，并且每个元素也需要加权。即
v
a
r
(
x
)
=
∑
w
i
(
x
i
−
u
x
′
)
2
var(x)=\sum w_i(x_i -u_x^{'})^2
var(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;
end
if (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;
return
end
window = fspecial('gaussian', 11, 1.5); %
K(1) = 0.01; % default settings
K(2) = 0.03; %
L = 255; %
end
if (nargin == 3)
if ((M < 11) || (N < 11))
ssim_index = -Inf;
ssim_map = -Inf;
return
end
window = 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;
end
else
ssim_index = -Inf;
ssim_map = -Inf;
return;
end
end
if (nargin == 4)
[H W] = size(window);
if ((H*W) < 4 || (H > M) || (W > N))
ssim_index = -Inf;
ssim_map = -Inf;
return
end
L = 255;
if (length(K) == 2)
if (K(1) < 0 || K(2) < 0)
ssim_index = -Inf;
ssim_map = -Inf;
return;
end
else
ssim_index = -Inf;
ssim_map = -Inf;
return;
end
end
if (nargin == 5)
[H W] = size(window);
if ((H*W) < 4 || (H > M) || (W > N))
ssim_index = -Inf;
ssim_map = -Inf;
return
end
if (length(K) == 2)
if (K(1) < 0 || K(2) < 0)
ssim_index = -Inf;
ssim_map = -Inf;
return;
end
else
ssim_index = -Inf;
ssim_map = -Inf;
return;
end
end
C1 = (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));
else
numerator1 = 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);
end
mssim = 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;
end
if 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);
end
else
% General case: Equation 12 from [1]
% Luminance term
if (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 term
sigmaxsigmay = [];
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 term
if (exponents(3) > 0)
num = sigmaxy + C(3);
if isempty(sigmaxsigmay)
sigmaxsigmay = sqrt(sigmax2.*sigmay2);
end
den = sigmaxsigmay + C(3);
ssimmap = ssimmap.*guardedDivideAndExponent(num,den,C(3),exponents(3));
end
end
ssimval = mean(ssimmap(:));
end
% -------------------------------------------------------------------------
function component = guardedDivideAndExponent(num, den, C, exponent)
if C > 0
component = num./den;
else
component = ones(size(num),'like',num);
isDenNonZero = (den ~= 0);
component(isDenNonZero) = num(isDenNonZero)./den(isDenNonZero);
end
if (exponent ~= 1)
component = component.^exponent;
end
end
function gaussFilt = getGaussianWeightingFilter(radius,N)
% Get 2D or 3D Gaussian weighting filter
filtRadius = 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
end
end
function [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'));
end
if ~isequal(size(A),size(ref))
error(message('images:validate:unequalSizeMatrices','A','REF'));
end
if (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:nargin
arg = 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) > 1
error(message('images:validate:ambiguousInputString', arg))
elseif numel(idx) == 1
if (i+1 > nargin)
error(message('images:validate:missingParameterValue'));
end
if idx == 1
dynmRange = varargin{i+1};
validateattributes(dynmRange,{'numeric'},{'positive', ...
'finite', 'real', 'nonempty','scalar'}, mfilename, ...
'DynamicRange',i);
dynmRange = double(dynmRange);
elseif idx == 2
C = varargin{i+1};
validateattributes(C,{'numeric'},{'nonnegative','finite', ...
'real','nonempty','vector', 'numel', 3}, mfilename, ...
'RegularizationConstants',i);
C = double(C);
elseif idx == 3
exponents = varargin{i+1};
validateattributes(exponents,{'numeric'},{'nonnegative', ...
'finite', 'real', 'nonempty','vector', 'numel', 3}, ...
mfilename,'Exponents',i);
exponents = double(exponents);
elseif idx == 4
radius = varargin{i+1};
validateattributes(radius,{'numeric'},{'positive','finite', ...
'real', 'nonempty','scalar'}, mfilename,'Radius',i);
radius = double(radius);
end
end
else
error(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];
end
end
对于代码核心部分的理解
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;
u x ′ u_x' ux′为加权均值,
v a r ( x ) = ∑ w i ( x i − u x ′ ) 2 = ∑ w i ( x i 2 − 2 u x ′ x i + u x ′ 2 ) = w i x i 2 − 2 u x ′ ∑ w i x i + ∑ w i u x ′ 2 = w i x i 2 + u x ′ ( u x ′ ∑ w i − 2 ∑ w i x i ) = ∑ w i x i 2 + u x ′ ( u x ′ − 2 u x ′ ) = ∑ w i x i 2 − u x ′ 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
注意 ∑ w i = 1 \sum w_i=1 ∑wi=1。
c o v ( x , y ) = ∑ w i ( x i − u x ′ ) ( y i − u y ′ ) = ∑ w i ( x i y i − x i u y ′ − u x ′ y i + u x ′ u y ′ ) = ∑ w i x i y i + u x ′ u y ′ − u y ′ ∑ w i x i − u x ′ ∑ w i y i = ∑ w i x i y i + u x ′ u y ′ − u y ′ u x ′ − u x ′ u y ′ = ∑ w i x i y i − u x ′ u y ′ \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′
438
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
