图像融合质量评价方法SSIM、PSNR、EN、MSE与NRMSE（一）

1 前言

最近开始学习和研究医疗图像融合这一方向，在做实验时因为缺少评估融合图像的指标以及代码非常苦恼。在经过一段时间的收集之后，得到了不少指标，所以希望可以记录一下。本文主要记录怎么使用这些指标，对于这些指标的公式和原理并不会过多介绍。后续更新图像融合质量评价方法AG、SF、STD、MI与NMI（二） 、图像融合质量评价方法MSSIM、MS-SSIM、FS、Qmi、Qabf与VIFF（三） 和图像融合质量评价方法FMI（四）

2 融合评价指标介绍

2.1 结构相似性 SSIM

SSIM是Structural Similarity的缩写，表示结构相似性，取值范围为[-1,1]，越接近1，代表相似度越高，融合质量越好。
关于详细的公式和原理解释可以看Image Quality Assessment: From Error Visibility to Structural Similarity。

2.2 峰值信噪比 PSNR

Peak signal to noise ration缩写PSNR，峰值信噪比，用于衡量图像有效信息与噪声之间的比率，能够反映图像是否失真。
$$PSNR=10\log \frac{z^2}{MSE}$$Z表示理想参考图像灰度最大值与最小值之差，通常为255。PSNR的值越大，表示融合图像的质量越好。

2.3 信息熵 Entropy

信息熵主要是度量图像包含信息量多少的一个客观评价指标,a表示灰度值，P(a)表示灰度概率分布:
$$H(A)=-\sum _a{P}_A(a)l\mathrm{o}\mathrm{g}{p}_A(a)$$信息熵越高表示融合图像的信息量越丰富，质量越好

2.4 均方误差 MSE

均方误差(Mean Square Error)反映的是变量间的差异程度，是一种基于像素误差的图像质量客观评价指标，用于衡量融合图像和理想参考图像之间的差异，MSE越小，表示融合图像质量越好。
均方根误差(Root Mean Square Error)是一个翻译空间细节信息的评价指标，公式如下：$$RMSE(I_F,I_W)=\sqrt{MSE(I_F,I_W)}$$

2.5 归一化均方根误差 NRMSE

归一化均方根误差（normalized root mean square error）就是将RMSE的值变成(0,1)之间。

3 代码实现

3.1 使用skimage库

为什么要一次性介绍这5个指标？因为这5个指标都只需要调用skimage库就可以使用了。skimage库是用于图像处理的 Python 包，使用原生的 NumPy 数组作为图像对象。
其中skimage.measure中包括以上所以指标的函数，详见Scikit image 中文开发手册。废话不多说，下面之间给出使用这些函数的例子：

import cv2
import numpy as np
import math

import skimage.measure 

# 参数0表示读入即为灰度图
img = cv2.imread('1.png',0)
img_fu = cv2.imread('2.png',0)

# 结构相似性
ssim = skimage.measure.compare_ssim(img, img_fu, data_range=255)

# 峰值信噪比
psnr = skimage.measure.compare_psnr(img, img_fu, data_range=255)

# 均方误差
mse = skimage.measure.compare_mse(img, img_fu)

# 均方根误差
rmse = math.sqrt(mse)

# 归一化均方根误差
nrmse = skimage.measure.compare_nrmse(img, img_fu, norm_type='Euclidean')

# 信息熵
entropy = skimage.measure.shannon_entropy(img_fu, base=2)

3.2 信息熵Entropy python实现

其实信息熵的实现并不难，根据上面的公式就可以写出（代码参考图像信息熵）：

import cv2
import numpy as np
import math

import skimage.measure 


def imageEn(image):
    tmp = []
    for i in range(256):
        tmp.append(0)
    val = 0
    k = 0
    res = 0
    img = np.array(image)
    for i in range(len(img)):
        for j in range(len(img[i])):
            val = img[i][j]
            tmp[val] = float(tmp[val] + 1)
            k = float(k + 1)
    for i in range(len(tmp)):
        tmp[i] = float(tmp[i] / k)
    for i in range(len(tmp)):
        if(tmp[i] == 0):
            res = res
        else:
            res = float(res - tmp[i] * (math.log(tmp[i]) / math.log(2.0)))

    return res

if __name__ == '__main__':
    image = cv2.imread('1.png',0)
    print(imageEn(image))
    print(skimage.measure.shannon_entropy(image, base=2))

可以看到对同一图片，上面所示代码和skimage库中的shannon_entropy()函数计算结果相差不大

3.3 峰值信噪比PSNR python实现

代码来源PSNR实现：

import cv2
import numpy as np
import math

import skimage.measure 

 
def psnr1(img1, img2):
   mse = np.mean((img1/1.0 - img2/1.0) ** 2 )
   if mse < 1.0e-10:
      return 100
   return 10 * math.log10(255.0**2/mse)
 
def psnr2(img1, img2):
   mse = np.mean( (img1/255. - img2/255.) ** 2 )
   if mse < 1.0e-10:
      return 100
   PIXEL_MAX = 1
   return 20 * math.log10(PIXEL_MAX / math.sqrt(mse))

if __name__ == '__main__':
    image = cv2.imread('1.png',0)
    img_fu = cv2.imread('2.png',0)
    print(psnr1(image,img_fu))
    print(psnr2(image,img_fu))
    print(skimage.measure.compare_psnr(image, img_fu, data_range=255))

对同一图片进行测试，得到结果如下：

3.4 结构相似性SSIM python实现

在研究使用python实现MS-SSIM（多层级结构相似性）时，考虑从SSIM代码上面进行改进，发现【超分辨率】基于结构相似性的图像质量评价SSIM及其python实现这篇文章的代码误差相对而言较小：

import numpy as np
from PIL import Image 
from scipy.signal import convolve2d
import cv2

import skimage.measure 

def matlab_style_gauss2D(shape=(3,3),sigma=0.5):
    """
    2D gaussian mask - should give the same result as MATLAB's
    fspecial('gaussian',[shape],[sigma])
    """
    m,n = [(ss-1.)/2. for ss in shape]
    y,x = np.ogrid[-m:m+1,-n:n+1]
    h = np.exp( -(x*x + y*y) / (2.*sigma*sigma) )
    h[ h < np.finfo(h.dtype).eps*h.max() ] = 0
    sumh = h.sum()
    if sumh != 0:
        h /= sumh
    return h

def filter2(x, kernel, mode='same'):
    return convolve2d(x, np.rot90(kernel, 2), mode=mode)

def compute_ssim(im1, im2, k1=0.01, k2=0.03, win_size=11, L=255):

    if not im1.shape == im2.shape:
        raise ValueError("Input Imagees must have the same dimensions")
    if len(im1.shape) > 2:
        raise ValueError("Please input the images with 1 channel")

    M, N = im1.shape
    C1 = (k1*L)**2
    C2 = (k2*L)**2
    window = matlab_style_gauss2D(shape=(win_size,win_size), sigma=1.5)
    window = window/np.sum(np.sum(window))

    if im1.dtype == np.uint8:
        im1 = np.double(im1)
    if im2.dtype == np.uint8:
        im2 = np.double(im2)

    mu1 = filter2(im1, window, 'valid')
    mu2 = filter2(im2, window, 'valid')
    mu1_sq = mu1 * mu1
    mu2_sq = mu2 * mu2
    mu1_mu2 = mu1 * mu2
    sigma1_sq = filter2(im1*im1, window, 'valid') - mu1_sq
    sigma2_sq = filter2(im2*im2, window, 'valid') - mu2_sq
    sigmal2 = filter2(im1*im2, window, 'valid') - mu1_mu2

    ssim_map = ((2*mu1_mu2+C1) * (2*sigmal2+C2)) / ((mu1_sq+mu2_sq+C1) * (sigma1_sq+sigma2_sq+C2))

    return np.mean(np.mean(ssim_map))


if __name__ == "__main__":
    im1 = cv2.imread('1.png',0)
    im2 = cv2.imread('2.png',0)

    print(compute_ssim(np.array(im1),np.array(im2)))
    ssim = skimage.measure.compare_ssim(im1, im2, data_range=255)
    print(ssim)

可以看到使用上面的代码与skimage库计算SSIM误差存在但是相对而言较小

4 总结

这5个指标的实现比较容易，但是融合图像评价方法还有很多，继续慢慢总结…