MSE与PSNR的计算方法

MSE: Mean Square Error，均方误差

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

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

MSE：
$$MSE=\frac{1}{N}\sum _i^N(x_i-y_i{)}^2$$
这里虽然将$X,Y$表达成一维的形式，但事实上其尺寸可任意。在进行图像质量评估时，MSE越小表明待评估的图像质量越好。

PSNR：
$$PNSR=10{\log }_{10}\frac{L^2}{MSE}$$
可以看出，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

理念

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

概念准备

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

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

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

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

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

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

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

公式理解

基本元素

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

均值：
$${\mu }_x=\frac{1}{N}\sum _i^N{x}_i\text{\hspace{1em}(1)}$$
标准差：
$${\sigma }_x=(\frac{1}{N-1}\sum _i^N(x_i-{\mu }_x{)}^2{)}^{1/2}\text{\hspace{1em}(2)}$$
假设我们有另外一个图像块，像素数目也为N，其中的元素记为$y_i$，那么两个图像块减去均值后的余弦相似度公式为：
$$\begin{array}{rl}s(x,y)& =\frac{{\sum }_i^N(x_i-{\mu }_x)(y_i-{\mu }_y)}{({\sum }_i^N(x_i-{\mu }_x{)}^2{)}^{1/2}({\sum }_i^N(y_i-{\mu }_y{)}^2{)}^{1/2}}\\ & =\frac{\frac{1}{N-1}{\sum }_i^N(x_i-{\mu }_x)(y_i-{\mu }_y)}{(\frac{1}{N-1}{\sum }_i^N(x_i-{\mu }_x{)}^2{)}^{1/2}(\frac{1}{N-1}{\sum }_i^N(y_i-{\mu }_y{)}^2{)}^{1/2}}\\ & =\frac{{\sigma }_{xy}}{{\sigma }_x{\sigma }_y}\end{array}$$
其中${\sigma }_{xy}$表示协方差：
$${\sigma }_{xy}=\frac{1}{N-1}\sum _i^N(x_i-{\mu }_x)(y_i-{\mu }_y)\text{\hspace{1em}(3)}$$

三原则

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

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

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

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

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

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

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

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

亮度的相似度：
$$l(x,y)=\frac{2{\mu }_x{\mu }_y}{{\mu }_x^2+{\mu }_y^2}\stackrel{stabilization}{=}\frac{2{\mu }_x{\mu }_y+C_1}{{\mu }_x^2+{\mu }_y^2+C_1}\text{\hspace{1em}(4)}$$
对比度的相似度：
$$c(x,y)=\frac{2{\sigma }_x{\sigma }_y}{{\sigma }_x^2+{\sigma }_y^2}\stackrel{stabilization}{=}\frac{2{\sigma }_x{\sigma }_y+C_2}{{\sigma }_x^2+{\sigma }_y^2+C_2}\text{\hspace{1em}(5)}$$
结构的相似度则直接借用了余弦相似度的概念，前面我们已经写过了，现在直接将其罗列如下：
$$s(x,y)=\frac{{\sigma }_{xy}}{{\sigma }_x{\sigma }_y}\stackrel{stabilization}{=}\frac{{\sigma }_{xy}+C_3}{{\sigma }_x{\sigma }_y+C_3}\text{\hspace{1em}(6)}$$
上面的$C_1,C_2,C_3$是常数，用来使计算更加稳定，防止分母出现过小的情况。一般来说，这几个常数相比较$\mu ,\sigma$等变量应当是一个小值。

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

最终公式

有了（4）（5）（6）三个公式后，作者将SSIM定义成了如下形式：
$$SSIM(x,y)=[L(x,y){]}^{\alpha }[c(x,y){]}^{\beta }[s(x,y){]}^{\gamma }$$
$\alpha ,\beta ,\gamma$用于调节三部分的权重。在实际计算时， 通常取$\alpha =\beta =\gamma =1$。再令$c_3=C_2/2$，然后将（4）（5）（6）代入上式，（5）的分子和（6）的分母可以消去，于是可得：
$$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)}\text{\hspace{1em}(7)}$$
其中
$$C_1=(K_{1}L{)}^2\text{\hspace{1em}(8)}$$

$$C_2=(K_{2}L{)}^2\text{\hspace{1em}(9)}$$

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

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

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

程序计算方法

计算公式推导

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

方差公式：
$$\begin{array}{rl}va{r}_x={\sigma }_x^2& =\frac{1}{N-1}\sum _i^N(x_i-{\mu }_x{)}^2\\ & =\frac{1}{N-1}\sum _i^N(x_i^2-2x_i{\mu }_x+{\mu }_x^2)\\ & =\frac{1}{N-1}(\sum _i^N{x}_i^2-2{\mu }_x\sum _i^N{x}_i+\sum _i^N{\mu }_x^2)\end{array}$$
根据公式（1）可知
$$\sum _i^N{x}_i=N\cdot {\mu }_x$$
由于${\mu }_x$是一个标量，所以
$$\sum _i^N{\mu }_x^2=N\cdot {\mu }_x^2$$
在此基础上继续上面的推导：
$$\begin{array}{rl}va{r}_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)\\ & =\frac{1}{N-1}(\sum _i^N{x}_i^2-2N\cdot {\mu }_x^2+N\cdot {\mu }_x^2)\\ & =\frac{1}{N-1}(\sum _i^N{x}_i^2-N\cdot {\mu }_x^2)\\ & =\frac{N}{N-1}(\frac{1}{N}\sum _i^N{x}_i^2-{\mu }_x^2)\\ & =\frac{N}{N-1}({\mu }_{x^2}-{\mu }_x^2)\end{array}\text{\hspace{1em}(11)}$$
此时$\mu$可以看做一个算子，它的功能是对其脚标变量求平均。

同理可以推导协方差公式，文章中用${\sigma }_{xy}$表示：
$$\begin{array}{rl}{\sigma }_{xy}& =\frac{1}{N-1}\sum _i^N(x_i-{\mu }_x)(y_i-{\mu }_y)\\ & =\frac{1}{N-1}\sum _i^N(x_i{y}_i-x_i{\mu }_y-y_i{\mu }_x+{\mu }_x{\mu }_y)\\ & =\frac{1}{N-1}(\sum _i^N{x}_i{y}_i-{\mu }_y\sum _i^N{x}_i-{\mu }_x\sum _i^N{y}_i+\sum _i^N{\mu }_x{\mu }_y)\\ & =\frac{1}{N-1}(\sum _i^N{x}_i{y}_i-N\cdot {\mu }_y{\mu }_x-N\cdot {\mu }_x{\mu }_y+N\cdot {\mu }_x{\mu }_y)\\ & =\frac{1}{N-1}(\sum _i^N{x}_i{y}_i-N\cdot {\mu }_x{\mu }_y)\\ & =\frac{N}{N-1}({\mu }_{xy}-{\mu }_x{\mu }_y)\end{array}\text{\hspace{1em}(12)}$$

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

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

进而有：
$${\mu }_x=\frac{1}{N}\sum x_i=\frac{1}{N}\sum {\bar{x}}_i\cdot L=L\cdot {\mu }_{\bar{x}}$$

$$\begin{array}{rl}{\sigma }_x& =(\frac{1}{N-1}\sum _i^N(x_i-{\mu }_x{)}^2{)}^{1/2}\\ & =(\frac{1}{N-1}\sum _i^N(L\cdot {\bar{x}}_i-L\cdot {\mu }_{\bar{x}}{)}^2{)}^{1/2}\\ & =L\cdot (\frac{1}{N-1}\sum _i^N{\bar{x}}_i-{\mu }_{\bar{x}}{)}^2{)}^{1/2}\\ & =L\cdot {\sigma }_{\bar{x}}\end{array}$$

$$\begin{array}{rl}{\sigma }_{xy}& =\frac{1}{N-1}\sum _i^N(x_i-{\mu }_x)(y_i-{\mu }_y)\\ & =\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}})\\ & =L^2\cdot {\sigma }_{\bar{x}\bar{y}}\end{array}$$

有了以上铺垫后，可以把SSIM公式改写为：
$$\begin{array}{rl}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)}\\ & =\frac{(2L^2{\mu }_{\bar{x}}{\mu }_{\bar{y}}+K_1^2{L}^2)(2L^2{\sigma }_{\bar{x}\bar{y}}+K_2^2{L}^2)}{(L^2{\mu }_{\bar{x}}^2+L^2{\mu }_{\bar{y}}^2+K_1^2{L}^2)(L^2{\sigma }_{\bar{x}}^2+L^2{\sigma }_{\bar{y}}^2+K_2^2{L}^2)}\\ & =\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{array}\text{\hspace{1em}(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

图像中越亮的地方表明损失越小，越暗的地方表明损失越大。

# -*- 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()