图像相似性评价指标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)]α[c(x,y)]β[l(x,y)]γSSIM(x,y)={[l(x,y)]}^{\alpha}{[c(x,y)]}^{\beta}{[l(x,y)]}^{\gamma}SSIM(x,y)=[l(x,y)]α[c(x,y)]β[l(x,y)]γ

l(x,y)=2μxμy+C1μx2+μy2+C1l(x,y)=\frac{2\mu_x\mu_y+C_1}{\mu_x^2+\mu_y^2+C_1}l(x,y)=μx2+μy2+C12μxμy+C1

c(x,y)=2σxσy+C2σx2+σy2+C2c(x,y)=\frac{2\sigma_x\sigma_y+C_2}{\sigma_x^2+\sigma_y^2+C_2}c(x,y)=σx2+σy2+C22σxσy+C2

s(x,y)=σxy+C3σxσy+C3s(x,y)=\frac{\sigma_{xy}+C_3}{\sigma_x\sigma_y+C_3}s(x,y)=σxσy+C3σxy+C3

上式中:

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

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

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

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

————μx,μy,σx,σy\mu_x,\mu_y,\sigma_x,\sigma_yμx,μy,σx,σy分别是x,y的均值和方差。
—————σxy\sigma_{xy}σxy:是x,y的协方差，为E[(X−E(X))(Y−E(Y))]E[(X-E(X))(Y-E(Y))]E[(X−E(X))(Y−E(Y))]

————C1,C2,C3C_1,C_2,C_3C1,C2,C3：都为常数，用以保证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)。计算时，一般会将参数设置为α=β=β=1,C3=C22\alpha=\beta=\beta=1,C_3=\frac{C_2}{2}α=β=β=1,C3=2C2,SSIM计算可简化为：

SSIM(x,y)=(2μxμy+C1)(2σxy+C2)(μx2+μy2+C1)(σx2+σy2+C2)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)}SSIM(x,y)=(μx2+μy2+C1)(σx2+σy2+C2)(2μxμy+C1)(2σxy+C2)

1.3.Python代码实现

import sys
import numpy
from scipy import signal
from scipy import ndimagedef 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 = 11sigma = 1.5window = fspecial_gauss(size, sigma)K1 = 0.01K2 = 0.03L = 255 #bitdepth of imageC1 = (K1*L)**2C2 = (K2*L)**2mu1 = signal.fftconvolve(window, img1, mode='valid')mu2 = signal.fftconvolve(window, img2, mode='valid')mu1_sq = mu1*mu1mu2_sq = mu2*mu2mu1_mu2 = mu1*mu2sigma1_sq = signal.fftconvolve(window, img1*img1, mode='valid') - mu1_sqsigma2_sq = signal.fftconvolve(window, img2*img2, mode='valid') - mu2_sqsigma12 = signal.fftconvolve(window, img1*img2, mode='valid') - mu1_mu2if 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 = 5weight = numpy.array([0.0448, 0.2856, 0.3001, 0.2363, 0.1333])downsample_filter = numpy.ones((2, 2))/4.0im1 = 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=1ωh∑i=0ω−1∑j=0h−1[I(i,j)−K(i,j)]2MSE=\frac{1}{\omega h}\sum_{i=0}^{\omega-1}\sum_{j=0}^{h-1}[I(i,j)-K(i,j)]^2MSE=ωh1∑i=0ω−1∑j=0h−1[I(i,j)−K(i,j)]2

PSNR=10log10(MAXI2MSE)=20log10(MAX1MSE)PSNR=10{log}_{10}\left (\frac{MAX_I^2}{MSE} \right )=20{log}_{10}\left (\frac{MAX_1}{\sqrt{MSE}} \right )PSNR=10log10(MSEMAXI2)=20log10(MSEMAX1)

MAXIMAX_IMAXI是表示图像点颜色的最大数值，如果每个采样点用 8 位表示，则为28=2552^8=25528=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 - Imse = np.mean(np.square(error))psnr = 10 * np.log10(255 * 255 / mse)return psnrpsnr_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.shapegauss = np.random.normal(0, 50, (row,col))gauss = gauss.reshape(row,col)noisy = image + gaussreturn noisyimg = 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