背景

某日突发奇想在知乎搜索了一波图像质量评价，发现如下一段话"图像质量评价领域的科研人员纷纷表示‘质量评价就是灌水的’；研究生也因被导师强迫做该方向而恐慌：‘图像质量评价新生，开学研一，导师不让做项目，我该怎么做才能够保证自己将来找到一份好工作？’，好心人士献计献策“依你的情况最好的选择是不要读研了，直接去工作”；一个IQA方向快被老板逼疯了的人表示‘我觉得(图像质量评价)没什么意义’；大众对质量评价实际意义的质疑声此起彼伏：‘我们为什么要评价图像的质量？或者说，我们评价图像的质量有什么实际应用？首先，现阶段，学术上评价的都是各种数据库和失真图像，已提出的方法成百上千，然而即使是面对自然图像Challenge数据库仍然倍感无力；其次，人类的评价标准是很复杂的，数据库隐藏的评价体系只代表那些被调查者的评价，并不一定契合整个人类的视觉一致性。那么，既没有一个完全正确的参考标准，又没有推而广之的实际应用，每年成千上万篇的评价图像质量核心期刊有何意义？’”(以上内容均来自知乎)

实际上，图像质量评价绝不是什么洪水野兽，也非“就是灌水”，就我个人的角度出发，我认为该领域更像是为一部分图像领域所服务的研究，更偏重学术应用，而非工程应用。本质上，图像质量评价(Image Quality Assessment, IQA)在做的，其实就是在探索人眼视觉评估背后潜在的规律，比如人们会认为怎样的图像是“视觉体验差的”，如何量化的描述人眼的视觉评估，这其实是挺有意思的一件事，唯一的缺陷就是难以运用于于工程应用。而机器学习、深度学习、强化学习等当前热门的方向，也是被誉为“在探索人的思维模式和思考方式的科学研究”，但真要依赖计算机科学揭秘人脑思考的潜在规律，这一天不是我们这一代人能够见证到的。当前的“学习”，仍然只能在通过数据集去训练一个非线性模型以最小化损失函数误差这一步停滞不前。

还是回到本文的主题，图像、质量以及评估，为什么要对图像的质量进行评估，其实很容易解释：如果没有一套权威、统一的量化指标，那么大家做实验，发论文的时候就都会说“本文提出的方法优于 the state-of-the-art method(最先进的方法)”，因为没有衡量标准，大家就天高任鸟飞，随意放飞自我了。为了避免这种情况，同时也为了保证科学研究的严谨性，图像质量评价也就随之而生了。所以在前面，我会评价该领域为“该领域更像是为一部分图像领域所服务的研究，更偏重学术应用，而非工程应用”。毕竟工程应用不在乎图像质量的评估，用户觉得好那就好，而学术研究则更偏重于量化的评估，指标好就是好，没道理我正确率比你高、图像质量评估结果比你好，还得不到认可吧？那么，图像质量评估IQA就是这么一回事。当然，有些时候IQA也有工程应用上的作用，比如制定行业标准、国家标准的时候（比如博主所在学院的院长，制定了一个在印刷行业非常权威的标准，所以有些时候IQA也能为非学术相关的应用提供帮助），这时候IQA就能派上用场，但是总归来说，工业的应用对于IQA的需求是很少的。

全参考与盲评估

IQA分为两种情况来讨论，第一种情况是我们有参考图像(Reference Image)，称为全参考(Full Reference)图像质量评估；另一种是无参考(Blind)图像质量评估。其实还有一类叫半参考图像，但是因为没有限定信息条件，这一类研究都挺水的，什么叫做半参考，半参考又是参考多少，这些都没有统一的标准，甚至有一些所谓的半参考质量评估用的信息比全参考还多还夸张，那还叫什么半参考，参考就是参考，无参考就是无参考，这能届于中间么？抄袭就是抄袭，没抄袭就是没抄袭，半抄袭是什么鬼，是这么个道理吧。所以半参考通常是不被认可的领域或者被归于全参考的领域。

对于全参考评估，很好理解，就是我有一张标准无失真图像，以此进行图像质量评估，而盲图像评估就是没有参照图像的图像质量评估。下面将介绍一些最常见的图像质量评估方法。

全参考：图像相似性——色差(Color Difference, CD)

单纯从向量数值的角度出发，如果要度量两个向量间的相似度，比较容易想到的就是欧几里得距离(Euclidean Distance)或者角度差。而色差就是图像中的欧式距离，其中最典型的代表为均方误差(Mean Square Error, MSE)，求解的是失真图像与参考图像在每一个像素点上的距离度量之和再求平均值，很好理解，如下式所示：

MSE=\frac{1}{m\times n}\sum^m_{i=1}\sum^n_{j=1}\left( x\left(i,j \right)- y\left(i,j \right)\right)^2

其中， $m,n$ 表示图像在长、宽上的像素数， $x,y$ 分别表示失真图像与参考图像。

MSE本质上计算的是sRGB颜色空间下的色差，存在的显著问题是人眼色彩感知不一致以及忽略了图像的空间信息。

而CIELAB颜色空间比sRGB颜色空间能更好地贴近人眼色彩视觉感知，而L*u*v*颜色空间进一步改进了这一点，是比较权威的均匀颜色空间(Uniform Color Space)，即色彩在其中任意方向上的定量变化与人眼色彩感知相一致。

实现sRGB颜色空间至L*u*v*的变换如下：

由于sRGB为设备相关颜色空间，而L*u*v*为设备无关颜色空间，因此需要经过一个转换至XYZ颜色空间（设备无关）的变换作为过渡，如下：

\left[ \begin{matrix} X\\ Y \\Z \end{matrix} \right]= \left[ \begin{matrix} 0.412456 & 0.357576 & 0.180438\\0.212673 & 0.715152 & 0.072175 \\0.019334 & 0.119192 & 0.950304 \end{matrix} \right] \left[ \begin{matrix} V'_R\\ V'_G \\V'_B \end{matrix} \right]

其中

V'_c=\left\{
\begin{aligned}
\frac{V_c}{12.92}\ \ \ \ \ \ \ \ \ \ \ \ \ \ & ,\ if \ \ V_c\le 0.04045 \\
(\frac{V_c+0.055}{1.055})^{2.4} & ,\ if \ \ V_c>0.04045
\end{aligned}
\right.

其中， $V'_c=\frac{C}{255}$ ， $C=R,G,B$ ，R、G、B为sRGB三通道的值。需要指出，XYZ颜色空间需要预先选择参考白以确定转换矩阵，此处我们默认选择D65光源下的参考白。

而L*u*v*的计算又是基于Lu'v'颜色空间的，需要注意的是两者不是一个东西，Lu'v'颜色空间的计算转换如下：

\left\{
\begin{aligned}
L= \left\{
\begin{aligned}
903.3Y\ \ \ \ \ \ \ \ \ \ \ \ \ & ,\ if \ \ Y\le 0.008856 \\
116(Y)^{1/3}-16 & ,\ if \ \ Y>0.008856
\end{aligned}
\right.\\
u'=\frac{4X}{X+15Y+3Z}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\
v'=\frac{9Y}{X+15Y+3Z}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\end{aligned}
\right.

而u*v*的计算方式如下：

\left\{
\begin{aligned}
L^*=L\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\
u^*=13L^* \times (u'-u'_0) \\
v^*=13L^* \times (v'-v'_0) \
\end{aligned}
\right.

其中 $u'_0,v'_0$ 为标准参考白的 $u',v'$ 值，一般我们采用D65光源下的参考白，结果为 $(u'_{D65},v'_{D65})=(0.1978,0.4451)$ 。（D65光源下参考白的XYZ三刺激值为 $(0.9505, 1, 1.088)$ ，可以自己通过公式验算）

而L*u*v*色差就是L*u*v*颜色空间下的欧几里得距离。

全参考：图像相似性——峰值信噪比(Peak Signal to Noise Ratio, PSNR)

峰值信噪比（PSNR），一种评价图像的客观标准，PSNR是“Peak Signal to Noise Ratio”的缩写。peak的意思是顶点。而ratio的意思是比率的。整个意思就是到达噪声比率的顶点信号，psnr一般是用于最大值信号和背景噪音之间的一个工程项目。通常在经过影像压缩之后，通常输出的影像都会在某种程度与原始影像不同。为了衡量经过处理后的影像品质，我们通常会参考PSNR值来衡量某个处理程序能否令人满意。它是原图像与被处理图像之间的均方误差相对于(2^n-1)^2的对数值(信号最大值的平方，n是每个采样值的比特数)，它的单位是dB。

由于图像、影像本质上是多维信号，所以可以通过PSNR计算失真程度，它跟MSE的功能相似（因为是基于MSE求解的，感觉就是有点耍流氓），如下式：

PSNR=10\times \log_{10}\big ( \frac{MAX^2}{MSE} \big )

其中，MAX为图像的动态范围最大值，如果以默认[0,255]的范围代入求解，会出现负值，而归一化后([0,1]动态范围的图像)的结果比较正常，所以在计算时需要预先处理，在matlab中就需要对图像做一个./255的处理。

全参考：图像相似性——图像结构相似性(Structural SIMilarity index, SSIM)

下图揭示了SSIM在考虑图像结构、图像局部空域信息表现上的优越性，也可以说暴露了色差和PSNR的问题所在。如下图所示，圆环上的图片的MSE与PSNR值均相同，但从主观视觉的角度出发去评估，我们认为最上方的图片明显是质量高于下方几张的，但从MSE和PSNR的量化值来看，圆环上的几张失真图像的图像质量是相同的，这显然并不符合视觉上的感知。这是因为人眼对于邻近区域内的变化都更敏感，而MSE和PSNR都没有考虑空域信息上的信息。

SSIM的计算方法如下：

SSIM(x,y)=\frac{(2\mu_x\mu_y+c_1)(2\sigma_{xy}+c_2)}{(\mu_x^2+\mu_y^2+c_1)(\sigma_x^2+\sigma_y^2+c_1)}

其中， $\mu_x,\mu_y$ 分别表示图像 $x$ 与图像 $y$ 的像素均值，而 $\sigma_x^2$ 与 $\sigma_y^2$ 分别表示图像的方差。而 $c_1=(0.01T)^2$ , $c_2=(0.03T)^2$ ，T表示图像像素值的动态范围。

但需要指出的是，SSIM并非在任何情况下都优于色差和峰值信噪比，特定情况下各有千秋，这么说更合适。

我们可以从几组测试图像来看看SSIM与MSE的不同。

可以发现，SSIM index map更准确地捕获到了人眼感知的误差，对于失真图像的边缘的描述显然是比MSE更准确的。下图说明了同样的问题。

再来看一组图像的失真在MSE与SSIM下的量化结果，下图以图1作为参考图像，其中图2图3的图像显然是比图像456的视觉质量更高的，但光从MSE的角度来看，23456的MSE基本相同，在该测试图像中的SSIM显然体现出了非常大的优越性。

盲参考：色彩丰富度——色彩丰富度(Color Colorfulness Index, CCI)

图像相似性是全参考图像中的典型情况，因为只有我们有参考图像，才能评估图像的失真情况。而盲参考图像质量评价就需要视具体应用而定了，以博主比较熟悉的图像色彩领域而言，色彩丰富度CCI是一种较常用的图像评估手段，用以描述图像的色彩信息的丰富情况，如下式：

其中， $\bar{S}$ 表示图像饱和度均值，而 $\sigma$ 表示饱和度的标准差。

不难理解，用图像的色彩饱和度和色彩饱和度的扰动情况表示色彩丰富程度，一方面，色彩可以用色调、饱和度、亮度三个维度来描述，那么通过饱和度均值及其偏差的等权相加可以很好的描述色彩的丰富情况，不过也有一定的缺陷，缺乏了对于二者权重的考虑。

代码实战部分

MSE

% Matlab 2017a
function MSE=MSE_compute(img1,img2)
    [m,n,k]=size(img1);
    [m2,n2,k2]=size(img2);
    if m==m2 && n==n2 && k==k2
        MSE=sqrt(sum(sum(sum((img1-img2).^2)))./(m*n));
    end
end

Luv CD

function luvcd=luvcolordifference(img1,img2)

    WHITEPOINT=[0.950456,1,1.088754];
    WHITEPOINTU=(4*WHITEPOINT(1))./(WHITEPOINT(1)+15*WHITEPOINT(2)+3*WHITEPOINT(3));
    WHITEPOINTV=(9*WHITEPOINT(1))./(WHITEPOINT(1)+15*WHITEPOINT(2)+3*WHITEPOINT(3));
    [m,n,k]=size(img1);
    img1=img1+0.0000001*ones(m,n,k);
    L=zeros(m,n);
    xyz=rgb2xyz(img1./255,'WhitePoint','d65');
    u=4*xyz(:,:,1)./(xyz(:,:,1)+15.*xyz(:,:,2)+3.*xyz(:,:,3));
    v=9*xyz(:,:,2)./(xyz(:,:,1)+15.*xyz(:,:,2)+3.*xyz(:,:,3));
    for i=1:m
        for j=1:n
            if xyz(i,j,2)<=0.008856
                L(i,j)=903.3*xyz(i,j,2);
            else
                L(i,j)=116*(xyz(i,j,2)).^(1/3)-16;
            end
        end
    end
    uu=13.*L.*(u-WHITEPOINTU.*ones(m,n)); 	 
    vv=13.*L.*(v-WHITEPOINTV.*ones(m,n));


    img2=img2+0.0000001*ones(m,n,k);
    L2=zeros(m,n);
    xyz2=rgb2xyz(img2./255,'WhitePoint','d65');
    u2=4*xyz2(:,:,1)./(xyz2(:,:,1)+15.*xyz2(:,:,2)+3.*xyz2(:,:,3));
    v2=9*xyz2(:,:,2)./(xyz2(:,:,1)+15.*xyz2(:,:,2)+3.*xyz2(:,:,3));
    for i1=1:m
        for j1=1:n
            if xyz2(i1,j1,2)<=0.008856
                L2(i1,j1)=903.3*xyz2(i1,j1,2);
            else
                L2(i1,j1)=116*(xyz2(i1,j1,2)).^(1/3)-16;
            end
        end
    end
    uu2=13.*L2.*(u2-WHITEPOINTU.*ones(m,n)); 	 
    vv2=13.*L2.*(v2-WHITEPOINTV.*ones(m,n));


	luvcd = mean(mean(sqrt((L-L2).^2+(uu-uu2).^2 +(vv-vv2).^2)));
end

PSNR

function [peaksnr, snr] = psnr(A, ref, peakval)
checkImages(A,ref);
if nargin < 3
    peakval = diff(getrangefromclass(A));
else
    checkPeakval(peakval, A);
    peakval = double(peakval);
end

if isempty(A) % If A is empty, ref must also be empty
    peaksnr = zeros(0, 'like', A);
    snr     = zeros(0, 'like', A);
    return;
end

err = immse(A,ref);

peaksnr = 10*log10(peakval.^2/err);

if nargout > 1
    if isinteger(ref)  
        ref = double(ref);
    end                      
    snr = 10*log10(mean(ref(:).^2)/err);
end

end

function checkImages(A, ref)

validImageTypes = {'uint8','uint16','int16','single','double'};

validateattributes(A,validImageTypes,{'nonsparse'},mfilename,'A',1);
validateattributes(ref,validImageTypes,{'nonsparse'},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

end

function checkPeakval(peakval, A)

validateattributes(peakval,{'numeric'},{'nonnan', 'real', ...
    'nonnegative','nonsparse','nonempty','scalar'}, mfilename, ...
    'PEAKVAL',3);
if isinteger(A) && (peakval > diff(getrangefromclass(A)))
    warning(message('images:psnr:peakvalTooLarge', 'A', 'REF'));
end

end

SSIM

function [ssimval, ssimmap] = ssim(varargin)
narginchk(2,10);

[A, ref, C, exponents, radius] = parse_inputs(varargin{:});

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
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

CCI

img=double(imread('testimage.png'));
[m,n,k]=size(img);
img_hsv=rgb2hsv(img);
S_average=mean(mean(img_hsv(:,:,2)));
S_standarddeviation=sqrt(sum(sum(1/(m*n).*(img_hsv(:,:,2)-S_average.*ones(m,n)).^2)));
CCI=S_average+S_standarddeviation

function [h, s, v] = rgb2hsv(varargin)
[r, g, b, isColorMap, isEmptyInput, isThreeChannel] = parseInputs(varargin{:});

if(~isEmptyInput)

    if(isThreeChannel)
        imageIn(:,:,1) = r;
        imageIn(:,:,2) = g;
        imageIn(:,:,3) = b;
    elseif(isColorMap)
        imageIn = reshape(varargin{1},size(varargin{1},1),1,size(varargin{1},2));
    else
        imageIn = r;        
    end

    h = images.internal.rgb2hsvmex(imageIn);

    if(nargout == 3)
        s = h(:,:,2);
        v = h(:,:,3);
        h = h(:,:,1);        
    elseif(isColorMap)
        h = reshape(h,size(h,1), size(h,3));
    end
    
else
    if(isThreeChannel)
        h = r;
        s = g;
        v = b;
    else
        h = r;
    end
end
end

function [r, g, b, isColorMap, isEmptyInput, isThreeChannel] = parseInputs(varargin)

isColorMap = 0;
isEmptyInput = 0;
isThreeChannel = 0;

if (nargin == 1)
    r = varargin{1};
    g = [];
    b = [];
    if (ndims(r)==3)
        if isempty(r)
            isEmptyInput = 1;
            return
        end
        if(size(r,3) ~= 3)
            error(message('MATLAB:rgb2hsv:invalidInputSizeRGB'));
        end

        validateattributes(r, {'uint8', 'uint16', 'double', 'single'}, {'real'}, mfilename, 'RGB', 1);

    elseif ismatrix(r) %M x 3 matrix for M colors.
        
        isColorMap = 1;
        if(size(r,2) ~=3)
            error(message('MATLAB:rgb2hsv:invalidSizeForColormap'));
        end

        validateattributes(r, {'double'}, {'real','nonempty','nonsparse'}, mfilename, 'M');

        if((any(r(:) < 0) || any(r(:) > 1)))
            error(message('MATLAB:rgb2hsv:badMapValues'));
        end
        
    else
        error(message('MATLAB:rgb2hsv:invalidInputSize'));
    end

elseif (nargin == 3)
    isThreeChannel = 1;
    r = varargin{1};
    g = varargin{2};
    b = varargin{3};

    if isempty(r) || isempty(g) || isempty(b)
        isEmptyInput = 1;
        return
    end  
    
    validateattributes(r, {'uint8', 'uint16', 'double', 'single'}, {'real', '2d'}, mfilename, 'R', 1);
    validateattributes(g, {'uint8', 'uint16', 'double', 'single'}, {'real', '2d'}, mfilename, 'G', 2);
    validateattributes(b, {'uint8', 'uint16', 'double', 'single'}, {'real', '2d'}, mfilename, 'B', 3);

    if ~isa(r, class(g)) || ~isa(g, class(b)) || ~isa(r, class(b)) % mixed type inputs.
        r = im2double(r);
        g = im2double(g);
        b = im2double(b);
    end
    
    if ~isequal(size(r),size(g),size(b))
        error(message('MATLAB:rgb2hsv:InputSizeMismatch'));
    end
    
else
    error(message('MATLAB:rgb2hsv:WrongInputNum'));
end
end

实施例

clc,clear,close all
img1=double(imread('img_original.png'));
img2=double(imread('img_2.png'));
img3=double(imread('img_3.png'));
img4=double(imread('img_4.png'));
figure
subplot(221)
imshow(uint8(img1))
title('参考图像1')
subplot(222)
imshow(uint8(img2))
title('失真图像2')
subplot(223)
imshow(uint8(img3))
title('失真图像3')
subplot(224)
imshow(uint8(img4))
title('失真图像4')

载入图像并显示

%% MSE
MSEofimg2=MSE_compute(img1,img2)
MSEofimg3=MSE_compute(img1,img3)
MSEofimg4=MSE_compute(img1,img3)

MSE计算结果

MSEofimg2 =

   12.0006


MSEofimg3 =

   11.0482


MSEofimg4 =

   11.0482

可以发现三幅图像的MSE度量差距不大。

PSNR的计算如下，

%% PSNR
PSNRofimg2=psnr(img1./255,img2./255)
PSNRofimg3=psnr(img1./255,img3./255)
PSNRofimg4=psnr(img1./255,img4./255)

结果如下：

PSNRofimg2 =

   26.5468


PSNRofimg3 =

   27.2650


PSNRofimg4 =

   27.6282

同样地，PSNR的差别也不大。

SSIM：

%% SSIM
SSIMofimg2=ssim(img1./255,img2./255)
SSIMofimg3=ssim(img1./255,img3./255)
SSIMofimg4=ssim(img1./255,img4./255)

SSIMofimg2 =

    0.9871


SSIMofimg3 =

    0.7427


SSIMofimg4 =

    0.7162

在SSIM的度量下，终于能够比较明显的体现出三幅图像在视觉质量上的差异，需要注意的是，SSIM的求解也需要图像的像素动态范围在[0,1]之间。

由于色差主要应用在彩色图像计算上，因此我们不用先前的爱因斯坦的灰度图像进行色差计算，重新载入一对被改变了色彩的彩色图像。

clc,clear,close all
img1=double(imread('testcolorimg1.png'));
img2=double(imread('testcolorimg2.png'));
figure,subplot(121),imshow(uint8(img1)),title('原图像')
subplot(122),imshow(uint8(img2)),title('失真图像')

载入图像

%% 图像相似度指标计算
LuvCDofimg=luvcolordifference(img1,img2)
SSIMofimg=ssim(img1./255,img2./255)
MSE=MSE_compute(img1,img2)

LuvCDofimg =

   10.6537


SSIMofimg =

    0.6408
    
MSE =

   25.6127

最后是颜色丰富度的计算

img=img1;
[m,n,k]=size(img);
img_hsv=rgb2hsv(img);
S_average=mean(mean(img_hsv(:,:,2)));
S_standarddeviation=sqrt(sum(sum(1/(m*n).*(img_hsv(:,:,2)-S_average.*ones(m,n)).^2)));
CCI=S_average+S_standarddeviation

img=img2;
[m,n,k]=size(img);
img_hsv=rgb2hsv(img);
S_average=mean(mean(img_hsv(:,:,2)));
S_standarddeviation=sqrt(sum(sum(1/(m*n).*(img_hsv(:,:,2)-S_average.*ones(m,n)).^2)));
CCI2=S_average+S_standarddeviation

两幅图像的色彩丰富度计算结果如下：

CCI =

    0.9007


CCI2 =

    0.9491

量化评价表明测试彩色图像2（底色为蓝色的枫叶图像）比测试彩色图像1的色彩更丰富。

写在最后

图像质量评估是在图像领域中用以作证所提出方法优越性的指标，对于绝大多数的研究领域都应当存在相应的评价指标。因此如果你发现某一个领域缺乏相关的研究，那你不妨可以花一些时间在这方面，如果你能够提出具有充分的可解释性的量化指标，那你不仅可以发一篇论文证明你的科研能力，更能为后来的科研工作者提供极大的便利。

此外，对于某一些基于数据集的评估指标虽然与图像没有直接关联，但它们在图像领域有时也非常重要，有精力的朋友可以自行了解一下，比如查准率、查全率、K-折交叉验证。

其中K-折交叉验证是在数据集规模较小的时候采用的一种评价手段，非常实用，如下图所示

图像质量评价研究新秀：基于概率模型、基于神经网络、基于低级特征等图像质量评估方法

基于概率模型的方法：

这类方法首先建立图像特征与图像质量之间的统计概率模型, 大多采用多变量高斯分布描述概率分布. 对待评价图像, 提取特征后根据概率模型计算最大后验概率的图像质量, 或根据与概率模型的匹配程度(如特征间的距离) 估计图像质量。复杂的模型将需要更多的数据量, 这类方法只有当数据量较大时才可能取得较好的效果。

基于神经网络的方法：

这类方法先提取一定的图像变换域或空间特征, 再基于已知质量数据训练一个神经网络回归分析模型, 由图像特征预测图像质量。不过要我评价概率模型和神经网络的图像质量评估的话，我想说，基于训练（基于概率模型和神经网络的方法）的方法本质上都是在耍流氓= =一方面没法对人眼视觉特性的潜在规律进行揭露，另一方面你拿一个黑箱模型出来有什么用嘛....所以下面的方法是我最看好的！

基于低级特征的方法：

这类方法通过一些图像低级特征的失真程度评估图像的质量，如通过纹理、梯度以及颜色。其中有一些基于低级特征的方法目前仍是最先进的（缺陷在于计算速度慢，不过计算速度慢对于图像质量评估不是什么致命的缺陷）。

图像质量评估公开数据集

图像质量评价的数据库很多，各种失真类型针对各种图像，但公认度最高的是LIVE, CSIQ,TID2008和TID2013这四个数据集。这些库都提供了每幅失真图像的主观评分值，也就是图像质量量化参考值(Ground Truth)。前两个数据集都是针对常见失真类型为主，即加性高斯白噪声、高斯模糊、JPEG压缩和JPEG2000压缩，而TID2013包含失真图像数量有3000幅，主观实验打分人数是917人，权威性当然是最高的，但由于失真类型数量高达25种，同时也是最难的。下面两个表分别是数据库大致情况，和分别包括的失真类型。