图像处理中常用的相似度评估指标

导读

有时候我们想要计算两张图片是否相似，而用来衡量两张图片相似度的算法也有很多，例如：RMSE、PSNR、SSIM、UQI、SIFT以及深度学习等。这篇文章主要介绍，RMSE、PSNR、SSIM、UQI这些指标的计算和应用，关于SIFT算法来计算图像的相似度在后面的文章中再进行介绍

影响这些算法的结果也有很多因素，图片的噪声、平移、缩放、旋转、裁剪、透视变换等，都会影响到算法的最终结果，所以我们需要根据不同的应用场景来选择使用不同的算法

MSE和RMSE

MSE(mean squared error)：图像像素值的平方误差
RMSE(root mean squared error)：图像像素值的平方根误差

• 缺点

1. 对比两张图片的size必须要完全一致
2. 对缩放、旋转、裁剪敏感

• 优点

1. 计算速度快
2. 算法简单

• 适应场景

1. 视频中对比前后两帧的差别
2. 计算图片的相似

• 代码实现

def mean_squared_error(img1,img2):    #判断两张图片的shape是否一致    assert img1.shape == img2.shape,"images have different shape "                                     "{} and {}".format(img1.shape,img2.shape)    #计算两张图像的MSE    img_mse = np.mean((img1.astype(np.float64) - img2.astype(np.float64))**2)    #将图像的MSE量化为相似度    img_sim = img_mse / (255**2)    return img_mse,1 - img_simdef root_mean_squared_error(img1,img2):    #判断两张图片的shape是否一致    assert img1.shape == img2.shape,"images have different shape "                                     "{} and {}".format(img1.shape,img2.shape)    #计算两张图像的MSE    img_rmse = np.sqrt(np.mean((img1.astype(np.float64) - img2.astype(np.float64))**2))    #将图像的MSE量化为相似度    img_sim = img_rmse / 255    return img_rmse,1 - img_sim

计算两张图片的MSE和RMSE

print(mean_squared_error(img1,img2))#(10826.879833697936, 0.8334966576901509)print(root_mean_squared_error(img1,img2))#(104.05229374549096, 0.5919517892333688)

PSNR

PSNR(Peak Signal Noise Ratio)也叫峰值信噪比：为了衡量处理后图像的品质，我们经常会使用到PSNR来衡量程序的处理结果是否令人满意。它是原图像与被处理图像之间的均方误差相对于 ( 2 n − 1 ) 2 (2^n-1)^2 (2n−1)2的对数值(信号最大值的平方，n是每个采样值的比特数)，计算公式如下：
P S N R = 10 ∗ l o g 10 ( 2 n − 1 ) 2 M S E PSNR = 10 * log_{10}frac{(2^n-1)^2}{MSE} PSNR=10∗log10​MSE(2n−1)2​

• 缺点
  PSNR的分数无法和人眼看到的品质完全一致

• 优点
  算法简单，应用广泛

• 适用场景
  用来评价画质客观量的测法，也可以用于图像相似的计算

• 代码实现

def peak_signal_noise_ratio(img1,img2):    #判断两张图片的size是否一致    assert img1.shape == img2.shape,"images hava different shape "                                     "{} and {}".format(img1.shape,img2.shape)    #获取数据类型所表示的最大值    MAX = np.iinfo(img1.dtype).max    #计算两张图片的MSE    mse,_ = mean_squared_error(img1,img2)    #计算两张图片的PSNR    psnr = 10 * np.log10(MAX**2 / mse)    return psnrprint(peak_signal_noise_ratio(img1,img2))#49.24985796979349

注：由于PSNR的取值范围在(0,inf)，PSNR的值越大表示图像越相似。如果想要将PSNR指标量化为相似度，可以基于测试图片上设置一个最大的PSNR值，进行取最大操作来换算成相似度

SSIM

SSIM(structural similarity index)结构相似性：是一种用来衡量两张图像相似程度的指标。当两张图像其中有一张为无失真图像，另一张为失真后的图像，两张图像的结构相似性可以看做是失真图像的品质质量指标。结构相似性相对于峰值信噪比而言，结构相似性指标在图像品质的衡量上更能符合人眼对图像品质的判断。

以 x x x和 y y y两个信号为例，结构相似性的计算公式如下：
S S I M ( x , y ) = ( l ( x , y ) ) α ( c ( x , y ) ) β ( s ( x , y ) ) γ l ( x , y ) = 2 μ x μ y + C 1 μ x 2 + μ y 2 + C 1 c ( x , y ) = 2 δ x δ y + C 2 δ x 2 + δ y 2 + C 2 s ( x , y ) = δ x y + C 3 δ x δ y + C 3 begin{aligned} SSIM(x,y) &= (l(x,y))^alpha(c(x,y))^beta(s(x,y))^gamma l(x,y) &= frac{2mu_xmu_y + C_1}{mu_x^2+mu_y^2+C_1} c(x,y) &= frac{2delta_xdelta_y + C_2}{delta_x^2+delta_y^2+C_2} s(x,y) &= frac{delta_{xy}+C_3}{delta_xdelta_y+C_3} end{aligned} SSIM(x,y)l(x,y)c(x,y)s(x,y)​=(l(x,y))α(c(x,y))β(s(x,y))γ=μx2​+μy2​+C1​2μx​μy​+C1​​=δx2​+δy2​+C2​2δx​δy​+C2​​=δx​δy​+C3​δxy​+C3​​​
上式中的 l ( x , y ) l(x,y) l(x,y)比较的是两个信号之间的亮度， c ( x , y ) c(x,y) c(x,y)比较的是两个信号的对比度， s ( x , y ) s(x,y) s(x,y)比较的是两个信号的结构，其中 α > 0 ， β > 0 ， γ > 0 alpha > 0，beta > 0，gamma > 0 α>0，β>0，γ>0主要是用来调整 l ( x , y ) 、 c ( x , y ) 、 s ( x , y ) l(x,y)、c(x,y)、s(x,y) l(x,y)、c(x,y)、s(x,y)的重要性， μ x 和 μ y mu_x和mu_y μx​和μy​表示的是信号的平均值， δ x 和 δ y delta_x和delta_y δx​和δy​表示的是信号的标准差，而 δ x y delta_{xy} δxy​为信号的协方差， C 1 、 C 2 、 C 3 C_1、C_2、C_3 C1​、C2​、C3​皆为常数，用来维持 l ( x , y ) 、 c ( x , y ) 、 s ( x , y ) l(x,y)、c(x,y)、s(x,y) l(x,y)、c(x,y)、s(x,y)的稳定，SSIM的值越大，表示两个信号之间的相似度越高

在计算图像结构相似性的时候，我们一般会将参数设为 α = β = γ = 1 以 及 C 3 = C 2 / 2 alpha=beta=gamma=1以及C_3=C_2/2 α=β=γ=1以及C3​=C2​/2所以 S S I M SSIM 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 ) SSIM(x,y) = frac{(2mu_xmu_y+C_1)(2delta_{xy}+C_2))}{(mu_x^2+mu_y^2+C_1)(delta_x^2+delta_y^2+C_2)} SSIM(x,y)=(μx2​+μy2​+C1​)(δx2​+δy2​+C2​)(2μx​μy​+C1​)(2δxy​+C2​))​
通常在计算图像的 S S I M SSIM SSIM的时候，我们都会采用移动的filter来进行计算

• 代码实现

from scipy import signalfrom enum import Enumdef filter2(img,fltr,mode='same'):    return signal.convolve2d(img, np.rot90(fltr,2), mode=mode)def _get_sums(GT,P,win,mode='same'):    mu1,mu2 = (filter2(GT,win,mode),filter2(P,win,mode))    return mu1*mu1, mu2*mu2, mu1*mu2def _get_sigmas(GT,P,win,mode='same',**kwargs):    if 'sums' in kwargs:        GT_sum_sq,P_sum_sq,GT_P_sum_mul = kwargs['sums']    else:        GT_sum_sq,P_sum_sq,GT_P_sum_mul = _get_sums(GT,P,win,mode)    return filter2(GT*GT,win,mode)  - GT_sum_sq,            filter2(P*P,win,mode)  - P_sum_sq,             filter2(GT*P,win,mode) - GT_P_sum_mulclass Filter(Enum):    UNIFORM = 0    GAUSSIAN = 1def fspecial(fltr,ws,**kwargs):    if fltr == Filter.UNIFORM:        return np.ones((ws,ws))/ ws**2    elif fltr == Filter.GAUSSIAN:        x, y = np.mgrid[-ws//2 + 1:ws//2 + 1, -ws//2 + 1:ws//2 + 1]        g = np.exp(-((x**2 + y**2)/(2.0*kwargs['sigma']**2)))        g[ g < np.finfo(g.dtype).eps*g.max() ] = 0        assert g.shape == (ws,ws)        den = g.sum()        if den !=0:            g/=den        return g    return Nonedef _ssim_single(GT, P, ws, C1, C2, fltr_specs, mode):    win = fspecial(**fltr_specs)    GT_sum_sq, P_sum_sq, GT_P_sum_mul = _get_sums(GT, P, win, mode)    sigmaGT_sq, sigmaP_sq, sigmaGT_P = _get_sigmas(GT, P, win, mode, sums=(GT_sum_sq, P_sum_sq, GT_P_sum_mul))    assert C1 > 0    assert C2 > 0    ssim_map = ((2 * GT_P_sum_mul + C1) * (2 * sigmaGT_P + C2)) / (                (GT_sum_sq + P_sum_sq + C1) * (sigmaGT_sq + sigmaP_sq + C2))    cs_map = (2 * sigmaGT_P + C2) / (sigmaGT_sq + sigmaP_sq + C2)    return np.mean(ssim_map), np.mean(cs_map)def ssim(GT,P,ws=11,K1=0.01,K2=0.03,MAX=None,fltr_specs=None,mode='valid'):    if MAX is None:        MAX = np.iinfo(GT.dtype).max    assert GT.shape == P.shape, "Supplied images have different sizes " +                                 str(GT.shape) + " and " + str(P.shape)    if fltr_specs is None:        fltr_specs=dict(fltr=Filter.UNIFORM,ws=ws)    C1 = (K1*MAX)**2    C2 = (K2*MAX)**2    ssims = []    css = []    for i in range(GT.shape[2]):        ssim,cs = _ssim_single(GT[:,:,i],P[:,:,i],ws,C1,C2,fltr_specs,mode)        ssims.append(ssim)        css.append(cs)    return np.mean(ssims),np.mean(css)img_sim,_ = ssim(img1,img2)#0.9459787655432684

UQI

UQI(Universal Quality Image Index)也叫图像通用质量指标：它能够很容易的计算并且应用到各种图像处理的应用中，主要结合三个因素来计算，相关性的损失、亮度失真、对比度失真。尽管UQI指标是从数学的角度来定义的不是利用人类的视觉系统，但是通过各种各样的图像失真表明UQI与主观质量测量有着惊人的一致性。

UQI指标的计算公式如下：
Q = 4 δ x y x ^ y ^ ( δ x 2 + δ y 2 ) [ ( x ^ ) 2 + ( y ^ ) 2 ] x ^ = 1 N ∑ i = 1 N x i y ^ = 1 N ∑ i = 1 N y i δ x 2 = 1 N − 1 ∑ i = 1 N ( x i − x ^ ) 2 δ y 2 = 1 N − 1 ∑ i = 1 N ( y i − y ^ ) 2 begin{aligned} Q &= frac{4delta_{xy}hat{x}hat{y}}{(delta^2_x+delta^2_y)[(hat{x})^2+(hat{y})^2]} hat{x} &= frac{1}{N}sum_{i=1}^{N}x_i hat{y} &= frac{1}{N}sum_{i=1}^{N}y_i delta_{x}^2 &= frac{1}{N - 1}sum_{i=1}^{N}(x_i-hat{x}) ^2 delta_{y}^2 &= frac{1}{N - 1}sum_{i=1}^{N}(y_i - hat{y})^2 end{aligned} Qx^y^​δx2​δy2​​=(δx2​+δy2​)[(x^)2+(y^​)2]4δxy​x^y^​​=N1​i=1∑N​xi​=N1​i=1∑N​yi​=N−11​i=1∑N​(xi​−x^)2=N−11​i=1∑N​(yi​−y^​)2​
Q的取值范围在[-1,1]之间

• 代码实现

import numpy as npfrom scipy.ndimage import uniform_filterdef _uqi_single(GT,P,ws):    N = ws**2    window = np.ones((ws,ws))    GT_sq = GT*GT    P_sq = P*P    GT_P = GT*P    GT_sum = uniform_filter(GT, ws)    P_sum =  uniform_filter(P, ws)    GT_sq_sum = uniform_filter(GT_sq, ws)    P_sq_sum = uniform_filter(P_sq, ws)    GT_P_sum = uniform_filter(GT_P, ws)    GT_P_sum_mul = GT_sum*P_sum    GT_P_sum_sq_sum_mul = GT_sum*GT_sum + P_sum*P_sum    numerator = 4*(N*GT_P_sum - GT_P_sum_mul)*GT_P_sum_mul    denominator1 = N*(GT_sq_sum + P_sq_sum) - GT_P_sum_sq_sum_mul    denominator = denominator1*GT_P_sum_sq_sum_mul    q_map = np.ones(denominator.shape)    index = np.logical_and((denominator1 == 0) , (GT_P_sum_sq_sum_mul != 0))    q_map[index] = 2*GT_P_sum_mul[index]/GT_P_sum_sq_sum_mul[index]    index = (denominator != 0)    q_map[index] = numerator[index]/denominator[index]    s = int(np.round(ws/2))    return np.mean(q_map[s:-s,s:-s])def uqi (GT,P,ws=8):    if len(GT.shape) == 2:        GT = GT[:, :, np.newaxis]        P = P[:, :, np.newaxis]    GT = GT.astype(np.float64)    P = P.astype(np.float64)    return np.mean([_uqi_single(GT[:,:,i],P[:,:,i],ws) for i in range(GT.shape[2])])

总结

除了上面的图像相似度评估指标之外，Sewar还提供了很多其他的图像质量评估指标算法，例如：MS-SSIM、ERGAS、SCC、RASE、SAM、D_lambda、D_S、QNR、VIF以及PSNR-B。

下面我们来对比一下一些常用指标在添加不同噪声的情况下，计算出来的相似度结果之间的差异

给图片添加噪声可以参考我的这篇文章对图像添加(高斯/椒盐/泊松/斑点)噪声

参考

1. https://ece.uwaterloo.ca/~z70wang/publications/uqi.html
2. https://github.com/andrewekhalel/sewar
3. https://ece.uwaterloo.ca/~z70wang/publications/quality_2c.pdf