图像质量评估指标：MSE，PSNR，SSIM

谁践踏了优雅 2022-11-21 03:49 1316阅读 0赞

### 文章目录 ###

*  MSE与PSNR的计算方法
 *  MSE与PSNR的问题
 *  SSIM的理念与计算方法
 *   *  理念
     *  概念准备
     *  公式理解
     *   *  基本元素
         *  三原则
         *  亮度、对比度、结构的相似度指标设计
         *  最终公式
     *  程序计算方法
     *   *  计算公式推导
         *  程序

# MSE与PSNR的计算方法 #

MSE: Mean Square Error，均方误差

PSNR: Peak Signal to Noise Ratio，峰值信噪比

假设现在有两个图像，名字分别是 X , Y X, Y X,Y，其中的元素分别记为 x i , y i x\_i, y\_i xi,yi，Array中包含的元素个数为N，那么MSE和PSNR的计算公式分别为：

**MSE：**  
 M S E = 1 N ∑ i N ( x i − y i ) 2 MSE = \\frac \{1\}\{N\} \\sum\_\{i\}^\{N\}(x\_i - y\_i)^2 MSE=N1i∑N(xi−yi)2  
这里虽然将 X , Y X, Y X,Y表达成一维的形式，但事实上其尺寸可任意。在进行图像质量评估时，MSE越小表明待评估的图像质量越好。

**PSNR：**  
 P N S R = 10 log ⁡ 10 L 2 M S E PNSR = 10\\log\_\{10\}\{\\frac \{L^2\}\{MSE\}\} PNSR=10log10MSEL2  
可以看出，PSNR事实上只是将MSE换成了信号处理中常用的db（分贝）表达形式，所以PSNR与MSE没有什么本质上的不同，只是数字上看起来更加有区分度一些（从数学上讲这算不得什么根本性的东西，但从人的辨识度来讲，这个优势还是挺明显的）。由于PSNR公式将MSE放在了分母上，所以在进行图像质量评估时，PSNR数值越大表明待评估的图像质量越好，这一点与MSE相反。PSNR公式中的L是一个常数，它表示图像数据类型的最大动态范围。比如对于`float`型的图像数据，其取值范围是\[0, 1\]，所以L=1。对于`uint8`类型的图像数据，其取值范围是\[0, 255\]，所以L=255。一般情况下L就这两种情况，至于HDR（High Dynamic Range，宽动态范围）的图像，则根据图像的位深另算。

# MSE与PSNR的问题 #

MSE与PSNR的问题是，在计算每个位置上的像素差异时，其结果仅与当前位置的两个像素值有关，与其它任何位置上的像素无关。这也就是说，这种计算差异的方式仅仅将图像看成了一个个孤立的像素点，而忽略了图像内容所包含的一些视觉特征，特别是图像的**局部结构信息**。而图像质量的好坏极大程度上是一个主观感受，其中结构信息对人主观感受的影响非常之大。

因此，采用这种方式计算出来的差异有时不能很好地反映图像质量。这就是SSIM（Structural Similarity，结构相似性）所希望要解决的问题。

# SSIM的理念与计算方法 #

参考文献：[Image Quality Assessment : From Error Visibility to Structural Similarity，Wang et. al., 2004][Image Quality Assessment _ From Error Visibility to Structural Similarity_Wang et. al._ 2004]

## 理念 ##

SSIM是为了解决上述提到的MSE与PSNR的问题而提出的，所以SSIM在理念上，与后面二者最重要最根本的一个不同是，SSIM计算两张图像在每个位置上的差异时，不是在该位置上从两张图中各取一个像素，而是各取了**一个区域的像素**。个人认为这是大道，而具体的计算方法属于小道，不过也是设计精良的小道。

## 概念准备 ##

在具体的计算方面，SSIM认为两个图像块之间的差异由以下几个部分组成：

*  亮度，对应于均值
 *  对比度，对应于方差或标准差
 *  结构，对应于余弦相似度

上述三组对应关系如何理解？下面我们做一个简单的解释。

图像的**亮度**是图像像素值大小的一种描述。对于单个像素值而言，像素值越大，就越接近白色，也就越亮，反之越暗；对于整个图像而言，所有像素的平均值越大，整张图像就越亮，反之越暗。

图像的**对比度**是图像像素值在整个动态范围内分布情况的一种描述。图像像素值在动态范围内分布地越宽广，对比度就越好；反之图像像素值分布地越紧凑，对比度就差。而这种分布范围的描述恰恰就对应于方差或标准差的概念。

图像的**结构**指的是像素间相对位置以及比例关系的描述，它与图像像素值的绝对大小反而没什么关系。比如考虑这样一个简单的图像，图像中只画了一个正方形，如果不考虑数值类型的截断问题的话，那么我们对这张图像乘以任意数值，得到结果仍然是一个结构一模一样的正方形。这正是余弦相似度所表达的概念，余弦相似度又称为余弦夹角，是对两个向量方向差异的一种度量，它不考虑向量的大小。图像结构就可以看作是图像在向量空间中的方向。

这三个维度的信息形成了一种互补的关系，在三者的共同约束下形成了对图像质量的描述。比如余弦相似度可以描述图像结构的差异程度，但是它却无法保证图像绝对数值的差异，比如向量（1, 2, 3），（2, 4, 6），（3, 6, 9）三者两两之间的余弦相似度均为1，但是它们的数值大小和方差却各不相同。反过来看，图像的均值和方差相同的情况下也无法保证图像的结构相同，因为即使我们将图像的像素位置完全打乱，其均值和方差也会保持完全不变，但是结构却发生了翻天覆地得变化。

## 公式理解 ##

### 基本元素 ###

假设一个图像块的像素数目为N，其中的元素记为 x i x\_i xi（向量的维度不影响公式的理解，并且可以与原文保持一致）。那么图像块的均值、标准差（无偏）分别为：

**均值：**  
 μ x = 1 N ∑ i N x i (1) \\mu\_\{x\} = \\frac\{1\}\{N\}\\sum\_\{i\}^\{N\}\{x\_i\} \\tag\{1\} μx=N1i∑Nxi(1)  
**标准差：**  
 σ x = ( 1 N − 1 ∑ i N ( x i − μ x ) 2 ) 1 / 2 (2) \\sigma\_\{x\} = \\bigg( \\frac\{1\}\{N-1\}\\sum\_\{i\}^\{N\}(x\_i-\\mu\_\{x\})^2 \\bigg)^\{1/2\} \\tag\{2\} σx=(N−11i∑N(xi−μx)2)1/2(2)  
假设我们有另外一个图像块，像素数目也为N，其中的元素记为 y i y\_i yi，那么两个图像块减去均值后的**余弦相似度**公式为：  
 s ( x , y ) = ∑ i N ( x i − μ x ) ( y i − μ y ) ( ∑ i N ( x i − μ x ) 2 ) 1 / 2 ( ∑ i N ( y i − μ y ) 2 ) 1 / 2 = 1 N − 1 ∑ i N ( x i − μ x ) ( y i − μ y ) ( 1 N − 1 ∑ i N ( x i − μ x ) 2 ) 1 / 2 ( 1 N − 1 ∑ i N ( y i − μ y ) 2 ) 1 / 2 = σ x y σ x σ y \\begin\{aligned\} s(x,y) &= \\frac \{\\sum\_\{i\}^\{N\}(x\_i-\\mu\_\{x\})(y\_i-\\mu\_\{y\})\} \{\\bigg(\\sum\_\{i\}^\{N\}(x\_i-\\mu\_\{x\})^2 \\bigg)^\{1/2\} \\bigg(\\sum\_\{i\}^\{N\}(y\_i-\\mu\_\{y\})^2 \\bigg)^\{1/2\}\} \\\\\[8ex\] &= \\frac \{\\frac\{1\}\{N-1\} \\sum\_\{i\}^\{N\}(x\_i-\\mu\_\{x\})(y\_i-\\mu\_\{y\})\} \{\\bigg(\\frac\{1\}\{N-1\}\\sum\_\{i\}^\{N\}(x\_i-\\mu\_\{x\})^2 \\bigg)^\{1/2\} \\bigg(\\frac\{1\}\{N-1\}\\sum\_\{i\}^\{N\}(y\_i-\\mu\_\{y\})^2 \\bigg)^\{1/2\}\} \\\\\[8ex\] &= \\frac \{\\sigma\_\{xy\}\}\{\\sigma\_\{x\} \\sigma\_\{y\}\} \\end\{aligned\} s(x,y)=(∑iN(xi−μx)2)1/2(∑iN(yi−μy)2)1/2∑iN(xi−μx)(yi−μy)=(N−11∑iN(xi−μx)2)1/2(N−11∑iN(yi−μy)2)1/2N−11∑iN(xi−μx)(yi−μy)=σxσyσxy  
其中 σ x y \\sigma\_\{xy\} σxy表示**协方差**：  
 σ x y = 1 N − 1 ∑ i N ( x i − μ x ) ( y i − μ y ) (3) \\sigma\_\{xy\} = \\frac\{1\}\{N-1\} \\sum\_\{i\}^\{N\}(x\_i-\\mu\_\{x\})(y\_i-\\mu\_\{y\}) \\tag\{3\} σxy=N−11i∑N(xi−μx)(yi−μy)(3)

### 三原则 ###

本文所讨论的图像质量评估均指有标签（Ground Truth）情况下的图像质量评估，此时待评估图像的质量使用它与标签的相似度进行衡量。

作者认为，相似度指标的设计均应满足以下三个原则：

1.  对称性（Symmetry）： S i m ( X , Y ) = S i m ( Y , X ) Sim(X,Y) = Sim(Y,X) Sim(X,Y)=Sim(Y,X)
2.  存在边界（Boundedness）： S i m ( X , Y ) ≤ 1 Sim(X,Y) \\le 1 Sim(X,Y)≤1
3.  最大值条件唯一（Unique maximum）：当且仅当 X = Y X=Y X=Y时， S i m ( X , Y ) = 1 Sim(X,Y)=1 Sim(X,Y)=1成立 。 X = Y X=Y X=Y的离散形式是：对于所有的 i = 0 , 1 , 2 , . . . , N i=0, 1, 2, ..., N i=0,1,2,...,N，均有 x i = y i x\_i = y\_i xi=yi

以上三个条件很容易理解。

*  对称性不必讲了，很显然应该如此。
 *  我们通常对`相似程度`的衡量，最高就是100%，这也就是边界定为1的原因。
 *  当相似度为100%时，应当，也只能是两张图完全一样的情况。

### 亮度、对比度、结构的相似度指标设计 ###

在文章中，亮度和对比度的相似度指标设计使用了相同的原则，均使用了我们高中时期学的一个很简单的不等式关系：  
 a 2 + b 2 ≥ 2 a b a^2 + b^2 \\ge 2ab a2\+b2≥2ab  
当且仅当 a = b a=b a=b时上式的等号关系成立。

将 a 2 + b 2 a^2 + b^2 a2\+b2放到分母上，并将亮度（均值）和对比度（标准差）代入其中，即可得到亮度和对比度的相似度指标。

亮度的相似度：  
 l ( x , y ) = 2 μ x μ y μ x 2 + μ y 2 = s t a b i l i z a t i o n 2 μ x μ y + C 1 μ x 2 + μ y 2 + C 1 (4) l(x,y) = \\frac \{2\\mu\_x\\mu\_y\}\{\\mu\_x^2 + \\mu\_y^2\} \\xlongequal\{stabilization\} \\frac \{2\\mu\_x\\mu\_y + C\_1\}\{\\mu\_x^2 + \\mu\_y^2 + C\_1\} \\tag\{4\} l(x,y)=μx2\+μy22μxμystabilizationμx2\+μy2\+C12μxμy\+C1(4)  
对比度的相似度：  
 c ( x , y ) = 2 σ x σ y σ x 2 + σ y 2 = s t a b i l i z a t i o n 2 σ x σ y + C 2 σ x 2 + σ y 2 + C 2 (5) c(x,y) = \\frac \{2\\sigma\_x\\sigma\_y\}\{\\sigma\_x^2 + \\sigma\_y^2\} \\xlongequal\{stabilization\} \\frac \{2\\sigma\_x\\sigma\_y + C\_2\}\{\\sigma\_x^2 + \\sigma\_y^2 + C\_2\} \\tag\{5\} c(x,y)=σx2\+σy22σxσystabilizationσx2\+σy2\+C22σxσy\+C2(5)  
结构的相似度则直接借用了余弦相似度的概念，前面我们已经写过了，现在直接将其罗列如下：  
 s ( x , y ) = σ x y σ x σ y = s t a b i l i z a t i o n σ x y + C 3 σ x σ y + C 3 (6) s(x,y) = \\frac \{\\sigma\_\{xy\}\}\{\\sigma\_\{x\} \\sigma\_\{y\}\} \\xlongequal\{stabilization\} \\frac \{\\sigma\_\{xy\} + C\_3\}\{\\sigma\_\{x\} \\sigma\_\{y\} + C\_3\} \\tag\{6\} s(x,y)=σxσyσxystabilizationσxσy\+C3σxy\+C3(6)  
上面的 C 1 , C 2 , C 3 C\_1, C\_2, C\_3 C1,C2,C3是常数，用来使计算更加稳定，防止分母出现过小的情况。一般来说，这几个常数相比较 μ , σ \\mu, \\sigma μ,σ等变量应当是一个小值。

很容易知道，（4）（5）（6）三个公式均满足前面提到的三个原则。

### 最终公式 ###

有了（4）（5）（6）三个公式后，作者将SSIM定义成了如下形式：  
 S S I M ( x , y ) = \[ L ( x , y ) \] α \[ c ( x , y ) \] β \[ s ( x , y ) \] γ SSIM(x,y) = \[L(x,y)\]^\\alpha \[c(x,y)\]^\\beta \[s(x,y)\]^\\gamma SSIM(x,y)=\[L(x,y)\]α\[c(x,y)\]β\[s(x,y)\]γ  
 α , β , γ \\alpha, \\beta, \\gamma α,β,γ用于调节三部分的权重。在实际计算时， 通常取 α = β = γ = 1 \\alpha = \\beta = \\gamma = 1 α=β=γ=1。再令 c 3 = C 2 / 2 c\_3 = C\_2 / 2 c3=C2/2，然后将（4）（5）（6）代入上式，（5）的分子和（6）的分母可以消去，于是可得：  
 S S I M ( x , y ) = ( 2 μ x μ y + C 1 ) ( 2 σ x y + C 2 ) ( μ x 2 + μ y 2 + C 1 ) ( σ x 2 + σ y 2 + C 2 ) (7) 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\_2)\} \\tag\{7\} SSIM(x,y)=(μx2\+μy2\+C1)(σx2\+σy2\+C2)(2μxμy\+C1)(2σxy\+C2)(7)  
其中  
 C 1 = ( K 1 L ) 2 (8) C\_1 = (K\_1L)^2 \\tag\{8\} C1=(K1L)2(8)

C 2 = ( K 2 L ) 2 (9) C\_2 = (K\_2L)^2 \\tag\{9\} C2=(K2L)2(9)

（7）就是最终用于计算SSIM的公式。

公式（8）（9）之所以这样写，是为了抵消掉图像数值类型的影响。根据文章，在实际计算时，取 K 1 = 0.01 , K 2 = 0.03 K\_1=0.01, K\_2=0.03 K1=0.01,K2=0.03。

**特别注意**：按照上述方法，我们计算得到了一个图像块的SSIM值，在实际计算时，一个图像块通常取成正方形，然后将该SSIM值赋给图像块的中心位置。然后滑动这个正方形块，可以得到其他所有位置的SSIM值，于是我们就得到了一个SSIM图。然而在实际中，我们对两张图像进行计算得到的通常是一个SSIM值，这个值在文章中称为mean SSIM（MSSIM），计算方法是对SSIM图取平均：  
 M S S I M ( X , Y ) = 1 M ∑ j M S S I M ( x j , y j ) (10) MSSIM(X,Y) = \\frac\{1\}\{M\} \\sum\_\{j\}^\{M\} SSIM(x\_j, y\_j) \\tag\{10\} MSSIM(X,Y)=M1j∑MSSIM(xj,yj)(10)

## 程序计算方法 ##

### 计算公式推导 ###

在真正写程序做计算时，为了让计算更加模块化，条理更清晰，部分公式的推导方式会跟原文不太一样，需要再推导一番。这里有一部分内容参考了skimage.metrics的计算方式。

**方差公式：**  
 v a r x = σ x 2 = 1 N − 1 ∑ i N ( x i − μ x ) 2 = 1 N − 1 ∑ i N ( x i 2 − 2 x i μ x + μ x 2 ) = 1 N − 1 ( ∑ i N x i 2 − 2 μ x ∑ i N x i + ∑ i N μ x 2 ) \\begin\{aligned\} var\_x = \\sigma\_x^2 &= \\frac \{1\}\{N-1\} \\sum\_\{i\}^\{N\} (x\_i - \\mu\_x)^2 \\\\\[4ex\] &= \\frac \{1\}\{N-1\} \\sum\_\{i\}^\{N\} (x\_i^2 - 2x\_i\\mu\_x + \\mu\_x^2) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} (\\sum\_\{i\}^\{N\}x\_i^2 - 2\\mu\_x \\sum\_\{i\}^\{N\} x\_i + \\sum\_\{i\}^\{N\} \\mu\_x^2) \\\\\[4ex\] \\end\{aligned\} varx=σx2=N−11i∑N(xi−μx)2=N−11i∑N(xi2−2xiμx\+μx2)=N−11(i∑Nxi2−2μxi∑Nxi\+i∑Nμx2)  
根据公式（1）可知  
 ∑ i N x i = N ⋅ μ x \\sum\_\{i\}^\{N\} x\_i = N \\cdot \\mu\_x i∑Nxi=N⋅μx  
由于 μ x \\mu\_x μx是一个标量，所以  
 ∑ i N μ x 2 = N ⋅ μ x 2 \\sum\_\{i\}^\{N\} \\mu\_x^2 = N \\cdot \\mu\_x^2 i∑Nμx2=N⋅μx2  
在此基础上继续上面的推导：  
 v a r x = 1 N − 1 ( ∑ i N x i 2 − 2 μ x ∑ i N x i + ∑ i N μ x 2 ) = 1 N − 1 ( ∑ i N x i 2 − 2 N ⋅ μ x 2 + N ⋅ μ x 2 ) = 1 N − 1 ( ∑ i N x i 2 − N ⋅ μ x 2 ) = N N − 1 ( 1 N ∑ i N x i 2 − μ x 2 ) = N N − 1 ( μ x 2 − μ x 2 ) (11) \\begin\{aligned\} var\_x &= \\frac \{1\}\{N-1\} (\\sum\_\{i\}^\{N\}x\_i^2 - 2\\mu\_x \\sum\_\{i\}^\{N\} x\_i + \\sum\_\{i\}^\{N\} \\mu\_x^2) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} (\\sum\_\{i\}^\{N\}x\_i^2 - 2N \\cdot \\mu\_x^2 + N \\cdot \\mu\_x^2) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} (\\sum\_\{i\}^\{N\}x\_i^2 - N \\cdot \\mu\_x^2) \\\\\[4ex\] &= \\frac \{N\}\{N-1\} (\\frac \{1\}\{N\} \\sum\_\{i\}^\{N\}x\_i^2 - \\mu\_x^2) \\\\\[4ex\] &= \\frac \{N\}\{N-1\} (\\mu\_\{x^2\} - \\mu\_x^2) \\end\{aligned\} \\tag\{11\} varx=N−11(i∑Nxi2−2μxi∑Nxi\+i∑Nμx2)=N−11(i∑Nxi2−2N⋅μx2\+N⋅μx2)=N−11(i∑Nxi2−N⋅μx2)=N−1N(N1i∑Nxi2−μx2)=N−1N(μx2−μx2)(11)  
此时 μ \\mu μ可以看做一个算子，它的功能是对其脚标变量求平均。

同理可以推导**协方差公式**，文章中用 σ x y \\sigma\_\{xy\} σxy表示：  
 σ x y = 1 N − 1 ∑ i N ( x i − μ x ) ( y i − μ y ) = 1 N − 1 ∑ i N ( x i y i − x i μ y − y i μ x + μ x μ y ) = 1 N − 1 ( ∑ i N x i y i − μ y ∑ i N x i − μ x ∑ i N y i + ∑ i N μ x μ y ) = 1 N − 1 ( ∑ i N x i y i − N ⋅ μ y μ x − N ⋅ μ x μ y + N ⋅ μ x μ y ) = 1 N − 1 ( ∑ i N x i y i − N ⋅ μ x μ y ) = N N − 1 ( μ x y − μ x μ y ) (12) \\begin\{aligned\} \\sigma\_\{xy\} &= \\frac \{1\}\{N-1\} \\sum\_\{i\}^\{N\} (x\_i - \\mu\_x) (y\_i - \\mu\_y) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} \\sum\_\{i\}^\{N\} (x\_iy\_i - x\_i\\mu\_y - y\_i\\mu\_x+ \\mu\_x \\mu\_y) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} (\\sum\_\{i\}^\{N\}x\_iy\_i - \\mu\_y\\sum\_\{i\}^\{N\}x\_i - \\mu\_x\\sum\_\{i\}^\{N\}y\_i+ \\sum\_\{i\}^\{N\}\\mu\_x \\mu\_y) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} (\\sum\_\{i\}^\{N\}x\_iy\_i - N \\cdot \\mu\_y\\mu\_x - N \\cdot \\mu\_x\\mu\_y+ N \\cdot \\mu\_x \\mu\_y) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} (\\sum\_\{i\}^\{N\}x\_iy\_i - N \\cdot \\mu\_x\\mu\_y) \\\\\[4ex\] &= \\frac \{N\}\{N-1\} (\\mu\_\{xy\} - \\mu\_x\\mu\_y) \\end\{aligned\} \\tag\{12\} σxy=N−11i∑N(xi−μx)(yi−μy)=N−11i∑N(xiyi−xiμy−yiμx\+μxμy)=N−11(i∑Nxiyi−μyi∑Nxi−μxi∑Nyi\+i∑Nμxμy)=N−11(i∑Nxiyi−N⋅μyμx−N⋅μxμy\+N⋅μxμy)=N−11(i∑Nxiyi−N⋅μxμy)=N−1N(μxy−μxμy)(12)

下面我们通过**归一化**，将公式中的L去掉。

如果记  x ˉ \\bar x xˉ 是  x x x 归一化后的结果，那么有：  
 x = L ⋅ x ˉ x = L \\cdot \\bar x x=L⋅xˉ

进而有：  
 μ x = 1 N ∑ x i = 1 N ∑ x ˉ i ⋅ L = L ⋅ μ x ˉ \\mu\_x = \\frac\{1\}\{N\} \\sum x\_i= \\frac\{1\}\{N\} \\sum \\bar x\_i \\cdot L = L \\cdot \\mu\_\{\\bar x\} \\\\\[6ex\] μx=N1∑xi=N1∑xˉi⋅L=L⋅μxˉ

σ x = ( 1 N − 1 ∑ i N ( x i − μ x ) 2 ) 1 / 2 = ( 1 N − 1 ∑ i N ( L ⋅ x ˉ i − L ⋅ μ x ˉ ) 2 ) 1 / 2 = L ⋅ ( 1 N − 1 ∑ i N x ˉ i − μ x ˉ ) 2 ) 1 / 2 = L ⋅ σ x ˉ \\begin\{aligned\} \\sigma\_x &= \\bigg( \\frac\{1\}\{N-1\}\\sum\_\{i\}^\{N\}(x\_i-\\mu\_\{x\})^2 \\bigg)^\{1/2\} \\\\\[4ex\] &= \\bigg( \\frac\{1\}\{N-1\}\\sum\_\{i\}^\{N\}(L \\cdot \\bar x\_i - L \\cdot \\mu\_\{\\bar x\})^2 \\bigg)^\{1/2\} \\\\\[4ex\] &= L \\cdot \\bigg( \\frac\{1\}\{N-1\}\\sum\_\{i\}^\{N\} \\bar x\_i - \\mu\_\{\\bar x\})^2 \\bigg)^\{1/2\} \\\\\[4ex\] &= L \\cdot \\sigma\_\{\\bar x\} \\\\\[6ex\] \\end\{aligned\} σx=(N−11i∑N(xi−μx)2)1/2=(N−11i∑N(L⋅xˉi−L⋅μxˉ)2)1/2=L⋅(N−11i∑Nxˉi−μxˉ)2)1/2=L⋅σxˉ

σ x y = 1 N − 1 ∑ i N ( x i − μ x ) ( y i − μ y ) = 1 N − 1 ∑ i N ( L ⋅ x ˉ i − L ⋅ μ x ˉ ) ( L ⋅ y ˉ i − L ⋅ μ y ˉ ) = L 2 ⋅ σ x ˉ y ˉ \\begin\{aligned\} \\sigma\_\{xy\} &= \\frac \{1\}\{N-1\} \\sum\_\{i\}^\{N\} (x\_i - \\mu\_x) (y\_i - \\mu\_y) \\\\\[4ex\] &= \\frac \{1\}\{N-1\} \\sum\_\{i\}^\{N\} (L \\cdot \\bar x\_i - L \\cdot \\mu\_\{\\bar x\}) (L \\cdot \\bar y\_i - L \\cdot \\mu\_\{\\bar y\}) \\\\\[4ex\] &= L^2 \\cdot \\sigma\_\{\\bar x \\bar y\} \\end\{aligned\} \\\\\[6ex\] σxy=N−11i∑N(xi−μx)(yi−μy)=N−11i∑N(L⋅xˉi−L⋅μxˉ)(L⋅yˉi−L⋅μyˉ)=L2⋅σxˉyˉ

有了以上铺垫后，可以把SSIM公式改写为：  
 S S I M ( x , y ) = ( 2 μ x μ y + C 1 ) ( 2 σ x y + C 2 ) ( μ x 2 + μ y 2 + C 1 ) ( σ x 2 + σ y 2 + C 2 ) = ( 2 L 2 μ x ˉ μ y ˉ + K 1 2 L 2 ) ( 2 L 2 σ x ˉ y ˉ + K 2 2 L 2 ) ( L 2 μ x ˉ 2 + L 2 μ y ˉ 2 + K 1 2 L 2 ) ( L 2 σ x ˉ 2 + L 2 σ y ˉ 2 + K 2 2 L 2 ) = ( 2 μ x ˉ μ y ˉ + K 1 2 ) ( 2 σ x ˉ y ˉ + K 2 2 ) ( μ x ˉ 2 + μ y ˉ 2 + K 1 2 ) ( σ x ˉ 2 + σ y ˉ 2 + K 2 2 ) (13) \\begin\{aligned\} 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\_2)\} \\\\\[4ex\] &= \\frac \{(2L^2\\mu\_\{\\bar x\}\\mu\_\{\\bar y\} + K\_1^2L^2) (2L^2\\sigma\_\{\\bar x \\bar y\} + K\_2^2L^2)\}\{(L^2\\mu\_\{\\bar x\}^2 + L^2\\mu\_\{\\bar y\}^2 + K\_1^2L^2) (L^2\\sigma\_\{\\bar x\}^2 + L^2\\sigma\_\{\\bar y\}^2 + K\_2^2L^2)\} \\\\\[4ex\] &= \\frac \{(2\\mu\_\{\\bar x\}\\mu\_\{\\bar y\} + K\_1^2) (2\\sigma\_\{\\bar x \\bar y\} + K\_2^2)\}\{(\\mu\_\{\\bar x\}^2 + \\mu\_\{\\bar y\}^2 + K\_1^2) (\\sigma\_\{\\bar x\}^2 +\\sigma\_\{\\bar y\}^2 + K\_2^2)\} \\end\{aligned\} \\tag\{13\} SSIM(x,y)=(μx2\+μy2\+C1)(σx2\+σy2\+C2)(2μxμy\+C1)(2σxy\+C2)=(L2μxˉ2\+L2μyˉ2\+K12L2)(L2σxˉ2\+L2σyˉ2\+K22L2)(2L2μxˉμyˉ\+K12L2)(2L2σxˉyˉ\+K22L2)=(μxˉ2\+μyˉ2\+K12)(σxˉ2\+σyˉ2\+K22)(2μxˉμyˉ\+K12)(2σxˉyˉ\+K22)(13)

可以看出，只要一开始的时候用L把图像数据归一化一下，后面就可以再也不用L了。另外我们也知道，对于图像数据而言，如果不考虑uint8类型的截断误差，那么不同数据类型（float，uint8）只是相差一个因子而已，而这个因子在SSIM公式中可以消掉，所以数据类型不影响SSIM的计算结果。然而因为uint8有截断误差，所以如果图像原本是float，转成uint8后，SSIM的结果将有少许差异。

### 程序 ###

下面的程序使用了两种方式计算MSE，PSNR，SSIM。一种是自己按照公式写代码，一种是直接调用了skimage的metrics模块，用于验证自己写的代码的正确性。

MSE和PSNR公式比较简单，没有那么多弯弯绕绕，所以两种方式计算出来的结果相同。

SSIM的结果就有差异了，差异来自两个原因。

1.  虽然都使用卷积的方式对图像块求滑动平均，但是我对图像边界的处理时填0，而skimage.metrics对边界的处理是镜像填充，这个无所谓哪个更好，我只是觉得填0是更加通用一些的操作；
2.  skimage.metrics的卷积存在很严重的截断误差（可以单步追进去看看它调用的uniform\_filter函数），导致的结果是，图像数据分别为float型和uint8型时，ssim的计算结果相差甚远，这显然不太合理。所以感觉要慎重使用skimage.metrics对SSIM的计算。

验证的例子是lena图像，一张是原图，一张是对图像先做两倍下采样，再做两倍上采样得到的结果，我们可以看看重采样会给图像带来多大的损失。

计算得到的mssim和ssim\_map如下：

*  一个是转成灰度图计算的结果：mssim = 0.9356
 *  一个是彩图计算的结果：mssim = 0.8965

图像中越亮的地方表明损失越小，越暗的地方表明损失越大。  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JieTE5ODc_size_16_color_FFFFFF_t_70_pic_center]  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JieTE5ODc_size_16_color_FFFFFF_t_70_pic_center 1]

# -*- coding: utf-8 -*-
    import cv2
    import numpy as np
    from skimage import metrics
    from scipy import signal
    
    
    def compute_mse(X, Y):
        """ compute mean square error of two images Parameters: ----------- X, Y: numpy array two images data Returns: -------- mse: float mean square error """
        X = np.float32(X)
        Y = np.float32(Y)
        mse = np.mean((X - Y) ** 2, dtype=np.float64)
        return mse
    
    
    def compute_psnr(X, Y, data_range):
        """ compute peak signal to noise ratio of two images Parameters: ----------- X, Y: numpy array two images data Returns: -------- psnr: float peak signal to noise ratio """
        mse = compute_mse(X, Y)
        psnr = 10 * np.log10((data_range ** 2) / mse)
        return psnr
    
    
    def compute_ssim(X, Y, win_size=7, data_range=None):
        """ compute structural similarity of two images Parameters: ----------- X, Y: numpy array two images data win_size: int window size of image patch for computing ssim of one single position data_range: int or float maximum dynamic range of image data type Returns: -------- mssim: float mean structural similarity ssim_map: numpy array (float) structural similarity map, same shape as input images """
        assert X.shape == Y.shape, "X, Y must have same shape"
        assert X.dtype == Y.dtype, "X, Y must have same dtype"
        assert win_size <= np.min(X.shape[0:2]), \
            "win_size should be <= shorter edge of image"
        assert win_size % 2 == 1, "win_size must be odd"
        if data_range is None:
            if 'float' in str(X.dtype):
                data_range = 1
            elif 'uint8' in str(X.dtype):
                data_range = 255
            else:
                raise ValueError(
                    'image dtype must be uint8 or float when data_range is None')
    
        X = np.squeeze(X)
        Y = np.squeeze(Y)
        if X.ndim == 2:
            mssim, ssim_map = _ssim_one_channel(X, Y, win_size, data_range)
        elif X.ndim == 3:
            ssim_map = np.zeros(X.shape)
            for i in range(X.shape[2]):
                _, ssim_map[:, :, i] = _ssim_one_channel(
                    X[:, :, i], Y[:, :, i], win_size, data_range)
            mssim = np.mean(ssim_map)
        else:
            raise ValueError("image dimension must be 2 or 3")
        return mssim, ssim_map
    
    
    def _ssim_one_channel(X, Y, win_size, data_range):
        """ compute structural similarity of two single channel images Parameters: ----------- X, Y: numpy array two images data win_size: int window size of image patch for computing ssim of one single position data_range: int or float maximum dynamic range of image data type Returns: -------- mssim: float mean structural similarity ssim_map: numpy array (float) structural similarity map, same shape as input images """
        X, Y = normalize(X, Y, data_range)
        C1 = 0.01 ** 2
        C2 = 0.03 ** 2
    
        num = win_size ** 2
        kernel = np.ones([win_size, win_size]) / num
        mean_map_x = convolve2d(X, kernel)
        mean_map_y = convolve2d(Y, kernel)
    
        mean_map_xx = convolve2d(X * X, kernel)
        mean_map_yy = convolve2d(Y * Y, kernel)
        mean_map_xy = convolve2d(X * Y, kernel)
    
        cov_norm = num / (num - 1)
        var_x = cov_norm * (mean_map_xx - mean_map_x ** 2)
        var_y = cov_norm * (mean_map_yy - mean_map_y ** 2)
        covar_xy = cov_norm * (mean_map_xy - mean_map_x * mean_map_y)
    
        A1 = 2 * mean_map_x * mean_map_y + C1
        A2 = 2 * covar_xy + C2
        B1 = mean_map_x ** 2 + mean_map_y ** 2 + C1
        B2 = var_x + var_y + C2
    
        ssim_map = (A1 * A2) / (B1 * B2)
        mssim = np.mean(ssim_map)
        return mssim, ssim_map
    
    
    def normalize(X, Y, data_range):
        """ convert dtype of two images to float64, and then normalize them by data_range Paramters: ---------- X, Y: numpy array two images data data_range: int or float maximum dynamic range of image data type Returns: -------- X, Y: numpy array two images """
        X = X.astype(np.float64) / data_range
        Y = Y.astype(np.float64) / data_range
        return X, Y
    
    
    def convolve2d(image, kernel):
        """ convolve single channel image and kernel Parameters: ----------- image: numpy array single channel image data kernel: numpy array kernel data Returns: -------- result: numpy array image data, same shape as input image """
        result = signal.convolve2d(image, kernel, mode='same', boundary='fill')
        return result
    
    
    def color_to_gray(image, normalization=False):
        """ convert color image to gray image Parameters: ----------- image: numpy array color image data normalization: bool whether to do nomalization Returns: -------- image: numpy array gray image data """
        image = cv2.cvtColor(image, code=cv2.COLOR_BGR2GRAY)
        if normalization and 'uint8' in str(image.dtype):
            image = image.astype(np.float64) / 255
        return image
    
    
    def resample(image, factor, interpolation):
        """ down sample and then upsample image by factor using a certain interpolation method Parameters: ----------- image: numpy array image data factor: sampling factor interpolation: enum interpolation type Returns: -------- image_resample: numpy array resampled image data """
        height, width = image.shape[0:2]
        downsample_size = (int(width / factor), int(height / factor))
        image_resample = cv2.resize(image,
                                    dsize=downsample_size,
                                    interpolation=interpolation)
        image_resample = cv2.resize(image_resample, dsize=(width, height),
                                    interpolation=interpolation)
        return image_resample
    
    
    if __name__ == '__main__':
        image_original = cv2.imread('lena512color.tiff')
        image_resample = resample(image_original, 2, cv2.INTER_CUBIC)
    
        # psnr
        psnr = compute_psnr(image_original, image_resample, 255)
        psnr_skimage = metrics.peak_signal_noise_ratio(image_original,
                                                       image_resample,
                                                       data_range=255)
        print(psnr)
        print(psnr_skimage)
    
        # ssim
        image_original = color_to_gray(image_original, normalization=False)
        image_resample = color_to_gray(image_resample, normalization=False)
    
        mssim, ssim_map = compute_ssim(image_original,
                                       image_resample)
        # set multichannel as False if using gray image, else True
        mssim_skimage = metrics.structural_similarity(image_original,
                                                      image_resample,
                                                      multichannel=False)
        print(mssim)
        print(mssim_skimage)
    
        cv2.imshow('original', image_original)
        cv2.imshow('resample', image_resample)
        cv2.imshow('ssim_map', ssim_map)
        cv2.waitKey(0)
        cv2.destroyAllWindows()

[Image Quality Assessment _ From Error Visibility to Structural Similarity_Wang et. al._ 2004]: http://atc.umh.es/gatcom/bin/oqam/Referencias/Wang_Bovik_Sheikh_Simoncelli_2004.pdf
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JieTE5ODc_size_16_color_FFFFFF_t_70_pic_center]: /images/20221120/59846f8f31084d0fa3873353bfd27d69.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JieTE5ODc_size_16_color_FFFFFF_t_70_pic_center 1]: https://img-blog.csdnimg.cn/20201030233048372.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JieTE5ODc=,size_16,color_FFFFFF,t_70#pic_center