图像相似性评价指标SSIM/PSNR

1.结构相似性指标SSIM

参考自维基百科SSIM

1.1介绍

结构相似性指标（英文：structural similarity index，SSIM index），是一种用以衡量两张数字图象相似性的指标。结构相似性在于衡量数字图像相邻像素的关联性，图像中相邻像素的关联性反映了实际场景中物体的结构信息。因此，在设计图像失真的衡量指标时，必须考虑结构性失真。SSIM指标于2004年提出1。但当图像出现位移、缩放、旋转（皆属于非结构性的失真）的情况无法有效的反映。

1.2计算

给定两个图像x,y，两者的结构相似性SSIM定义为:

$SSIM(x,y)={[l(x,y)]}^{\alpha }{[c(x,y)]}^{\beta }{[l(x,y)]}^{\gamma }$

$l(x,y)=\frac{2{\mu }_x{\mu }_y+C_1}{{\mu }_x^2+{\mu }_y^2+C_1}$

$c(x,y)=\frac{2{\sigma }_x{\sigma }_y+C_2}{{\sigma }_x^2+{\sigma }_y^2+C_2}$

$s(x,y)=\frac{{\sigma }_{xy}+C_3}{{\sigma }_x{\sigma }_y+C_3}$

上式中:

———— l(x,y)：比较两个图像的亮度，像素均值

———— c(x,y)：比较两个图像的对比度，像素方差

———— s(x,y)：比较两个图像的结构，协方差，像素间关系。

————$\alpha ,\beta ,\gamma$:大于零，调整l(x,y),c(x,y),s(x,y)相对重要性的参数。

————${\mu }_x,{\mu }_y,{\sigma }_x,{\sigma }_y$分别是x,y的均值和方差。
—————${\sigma }_{xy}$:是x,y的协方差，为$E[(X-E(X))(Y-E(Y))]$

————$C_1,C_2,C_3$：都为常数，用以保证l(x,y),c(x,y),s(x,y)的稳定

SSIM值越大代表相似度越高。

特点：

• 1）结构相似性指标是对称的,即SSIM(x,y)=SSIM(y,x)。
• 2）结构相似性指标的范围为0到1。

实际使用中，会使用滑动窗口，在NXN的小区块中分别计算视窗内的结构相似性指标，最后将所有的局部结构相似性指标值求平均，即为两张图像的SSIM值，也被称作MSSIM(Mean SSIm)。计算时，一般会将参数设置为$\alpha =\beta =\beta =1,C_3=\frac{C_2}{2}$,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_2)}$

1.3.Python代码实现

import sys
import numpy
from scipy import signal
from scipy import ndimage


def fspecial_gauss(size, sigma):
    x, y = numpy.mgrid[-size//2 + 1:size//2 + 1, -size//2 + 1:size//2 + 1]
    g = numpy.exp(-((x**2 + y**2)/(2.0*sigma**2)))
    return g/g.sum()


def ssim(img1, img2, cs_map=False):
    img1 = img1.astype(numpy.float64)
    img2 = img2.astype(numpy.float64)
    size = 11
    sigma = 1.5
    window = fspecial_gauss(size, sigma)
    K1 = 0.01
    K2 = 0.03
    L = 255 #bitdepth of image
    C1 = (K1*L)**2
    C2 = (K2*L)**2
    mu1 = signal.fftconvolve(window, img1, mode='valid')
    mu2 = signal.fftconvolve(window, img2, mode='valid')
    mu1_sq = mu1*mu1
    mu2_sq = mu2*mu2
    mu1_mu2 = mu1*mu2
    sigma1_sq = signal.fftconvolve(window, img1*img1, mode='valid') - mu1_sq
    sigma2_sq = signal.fftconvolve(window, img2*img2, mode='valid') - mu2_sq
    sigma12 = signal.fftconvolve(window, img1*img2, mode='valid') - mu1_mu2
    if cs_map:
        return (((2*mu1_mu2 + C1)*(2*sigma12 + C2))/((mu1_sq + mu2_sq + C1)*
                    (sigma1_sq + sigma2_sq + C2)), 
                (2.0*sigma12 + C2)/(sigma1_sq + sigma2_sq + C2))
    else:
        return ((2*mu1_mu2 + C1)*(2*sigma12 + C2))/((mu1_sq + mu2_sq + C1)*
                    (sigma1_sq + sigma2_sq + C2))

def mssim(img1, img2):
    """
    refer to https://github.com/mubeta06/python/tree/master/signal_processing/sp
    """
    level = 5
    weight = numpy.array([0.0448, 0.2856, 0.3001, 0.2363, 0.1333])
    downsample_filter = numpy.ones((2, 2))/4.0
    im1 = img1.astype(numpy.float64)
    im2 = img2.astype(numpy.float64)
    mssim = numpy.array([])
    mcs = numpy.array([])
    for l in range(level):
        ssim_map, cs_map = ssim(im1, im2, cs_map=True)
        mssim = numpy.append(mssim, ssim_map.mean())
        mcs = numpy.append(mcs, cs_map.mean())
        filtered_im1 = ndimage.filters.convolve(im1, downsample_filter, 
                                                mode='reflect')
        filtered_im2 = ndimage.filters.convolve(im2, downsample_filter, 
                                                mode='reflect')
        im1 = filtered_im1[::2, ::2]
        im2 = filtered_im2[::2, ::2]
    return (numpy.prod(mcs[0:level-1]**weight[0:level-1])*
                    (mssim[level-1]**weight[level-1]))

mssim_val = mssim(img, noise_img)
ssim_val = ssim(img, noise_img)
print(f"mssim_val: {mssim_val}, ssim_val: {ssim_val.mean()}")

2.峰值信噪比PSNR

2.1介绍

峰值信噪比(Peak Signal to Noise Ratio, PSNR)，表示的是信号的最大功率与噪声功率的比值。峰值信噪比越高，表示噪声影响越小；峰值信噪比越低，表示噪声影响越大。3单位是分贝dB，大于30dB，压缩前后图像质量肉眼很难看出区别。

2.2计算方式

通常噪声常使用均方误差（Mean Square Error,MSE）来描述。两个宽高为w,h的灰度图I\K，I为无噪声图像，K为I的带噪声近似，则：

$MSE=\frac{1}{\omega h}{\sum }_{i=0}^{\omega -1}{\sum }_{j=0}^{h-1}[I(i,j)-K(i,j){]}^2$

$PSNR=10{log}_{10}\left(\frac{MA{X}_I^2}{MSE}\right)=20{log}_{10}\left(\frac{MA{X}_1}{\sqrt{MSE}}\right)$

$MA{X}_I$是表示图像点颜色的最大数值，如果每个采样点用 8 位表示，则为$2^8=255$

若为彩色图像，通常有三种方法来计算：

• 1)分别计算 RGB 三个通道的 PSNR，然后取平均值。
• 2）MSE使RGB三通道MSE的平均
• 3)图片转化为YCbCr格式，然后只计算Y分量也就是亮度分量的 PSNR

针对超光谱图像，我们需要针对不同波段分别计算 PSNR，然后取平均值，这个指标称为 MPSNR。

2.3Python实现

import numpy as np
def get_psnr(I, K):
    error = K - I
    mse = np.mean(np.square(error))
    psnr = 10 * np.log10(255 * 255 / mse)
    return psnr

psnr_value = get_psnr(img, noise_img)
print(f"psnr_value： {psnr_value}")

补充：

上述计算使用的是添加了高斯噪声的图像和原图。

import cv2
import numpy as np
file_path = "img.png"

def get_gauss_noise_image(image):
        row,col= image.shape
        gauss = np.random.normal(0, 50, (row,col))
        gauss = gauss.reshape(row,col)
        noisy = image + gauss
        return noisy

img = cv2.imread(file_path, cv2.IMREAD_GRAYSCALE)
cv2.imwrite("gray_text.png", img)
noise_img = get_gauss_noise_image(img)
cv2.imwrite("gray_noise_text.png", noise_img)
# pnsr=12

参考：
[1]:https://www.cns.nyu.edu/pub/lcv/wang03-preprint.pdf
[2]:https://zh.wikipedia.org/wiki/%E7%B5%90%E6%A7%8B%E7%9B%B8%E4%BC%BC%E6%80%A7#cite_note-SSIM-1
[3]:https://zh.wikipedia.org/wiki/%E5%B3%B0%E5%80%BC%E4%BF%A1%E5%99%AA%E6%AF%94