结构相似形特征(SSIM)原理简介及python实现

结构相似形特征是图像全参考评价(FR-IQA)中经典的一个方法，由Zhou Wang等人在2004年发表的论文《Image Quality Assessment: From Error Visibility to Structural Similarity》中提出。作为一种全参考方法，需要同时利用原始图像(未失真)和失真图像。

一、结构相似性（structural similarity）

      自然图像具有极高的结构性，表现在图像的像素间存在着很强的相关性，尤其是在空间相似的情况下。这些相关性在视觉场景中携带着关于物体结构的重要信息。我们假设人类视觉系统（HSV）主要从可视区域内获取结构信息。所以通过探测结构信息是否改变来感知图像失真的近似信息。

    大多数的基于误差敏感度（error sensitivity）的质量评估方法(如MSE,PSNR)使用线性变换来分解图像信号，这不会涉及到相关性。我们要讨论的SSIM就是要找到更加直接的方法来比较失真图像和参考图像的结构。

二、SSIM指数

    物体表面的亮度信息与照度和反射系数有关，且场景中的物体的结构与照度是独立的，反射系数与物体有关。我们可以通过分离照度对物体的影响来探索一张图像中的结构信息。这里，把与物体结构相关的亮度和对比度作为图像中结构信息的定义。因为一个场景中的亮度和对比度总是在变化的，所以我们可以通过分别对局部的处理来得到更精确的结果。

                  

      由SSIM测量系统可得相似度的测量可由三种对比模块组成，分别为：亮度，对比度，结构。接下来我们将会对这三模块函数进行定义。

      首先，对于离散信号，我们以平均灰度来作为亮度测量的估计：

                                                       (1)

       亮度对比函数l（x,y）是关于的函数。

       然后，由测量系统知道要把平均灰度值从信号中去除，对于离散信号，可使用标准差来做对比度估量值。

                                        （2）

       对比度对比函数c(x,y)就是的函数。

       接下来，信号被自己的标准差相除，结构对比函数就被定义成和的函数。

       最后，三个对比模块组合成一个完整的相似测量函数：

                                    （3）

       

       S（x,y）应该满足以下三个条件：

       （1） 对称性：；

       （2） 有界性：；

       （3） 最大值唯一性：当且仅当x=y时，S（x,y）=1 。

      

       现在，我们定义三个对比函数。

       亮度对比函数：

                                                  （4）

       常数是为了避免接近0时造成系统的不稳定。

       特别的，我们选择，L为图像灰度级数，对于8-bit灰度图像，L=255，。公式（4）满足上述三个条件。

       对比度对比函数：

                                               （5）

       常数，且。公式（5）依然满足上述三个条件。

       结构对比函数：

                                                     （6）

       其中

                                         （7）

      最后把三个函数组合起来，得到SSIM指数函数：

                                   （8）

      这里，用来调整三个模块间的重要性。

      为了得到简化形式，设，得到：

                                 （9）

所以结构相似度指数从图像组成的角度将结构信息定义为独立于亮度、对比度的反映场景中物体结构的属性，并将失真建模为亮度、对比度和结构三个不同因素的组合。用均值作为亮度的估计，标准差作为对比度的估计，协方差作为结构相似程度的度量。 

而在实际应用中，一般采用高斯函数计算图像的均值、方差以及协方差，而不是采用遍历像素点的方式，以换来更高的效率。

 

三、SSIM指数应用于图像质量评估

   在图像质量评估之中，局部求SSIM指数的效果要好于全局。

1. 第一，图像的统计特征通常在空间中分布不均；
2. 第二，图像的失真情况在空间中也是变化的；
3. 第三，在正常视距内，人们只能将视线聚焦在图像的一个区域内，所以局部处理更符合人类视觉系统的特点；
4. 第四，局部质量检测能得到图片空间质量变化的映射矩阵，结果可服务到其他应用中。

     所以，在上述公式中，都加入了一个8*8的方形窗，并且逐像素的遍历整幅图片。每一步计算，和SSIM都是基于窗口内像素的，最终得到一个SSIM指数映射矩阵，由局部SSIM指数组成。然而，简单的加窗会使映射矩阵出现不良的“分块”效应。为解决这问题，我们使用11*11的对称高斯加权函数作为加权窗口，标准差为1.5，且

                                                 （10）

 

则的估计值表示为：

                                               （11）

                                  （12）

                               （13）

       应用这种加窗方法，映射矩阵就可展现出局部各向同性的性质。

       在这里，经过一些实验总结，我们把K1设为0.01，K2设为0.03，然后用平均SSIM指数作为整幅图像的估计质量评价：

                    （14）

其中X,Y为图像，为局部SSIM指数在映射中的位置，MN为局部窗口的数量。

 

四、SSIM算法实现：

1、TensorFlow自带函数：

如PSNR一样，SSIM这种常用计算函数也被tensorflow收编了，我们只需在tf中调用ssim就可以了：

tf.image.ssim(x, y, 255)

源代码如下：

def ssim(img1, img2, max_val):
  """Computes SSIM index between img1 and img2.

  This function is based on the standard SSIM implementation from:
  Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. (2004). Image
  quality assessment: from error visibility to structural similarity. IEEE
  transactions on image processing.

  Note: The true SSIM is only defined on grayscale.  This function does not
  perform any colorspace transform.  (If input is already YUV, then it will
  compute YUV SSIM average.)

  Details:
    - 11x11 Gaussian filter of width 1.5 is used.
    - k1 = 0.01, k2 = 0.03 as in the original paper.

  The image sizes must be at least 11x11 because of the filter size.

  Example:
  # Read images from file.
      im1 = tf.decode_png('path/to/im1.png')
      im2 = tf.decode_png('path/to/im2.png')
      # Compute SSIM over tf.uint8 Tensors.
      ssim1 = tf.image.ssim(im1, im2, max_val=255)

      # Compute SSIM over tf.float32 Tensors.
      im1 = tf.image.convert_image_dtype(im1, tf.float32)
      im2 = tf.image.convert_image_dtype(im2, tf.float32)
      ssim2 = tf.image.ssim(im1, im2, max_val=1.0)
      # ssim1 and ssim2 both have type tf.float32 and are almost equal.
    img1: First image batch.
    img2: Second image batch.
    max_val: The dynamic range of the images (i.e., the difference between the
      maximum the and minimum allowed values).

  Returns:
    A tensor containing an SSIM value for each image in batch.  Returned SSIM
    values are in range (-1, 1], when pixel values are non-negative. Returns
    a tensor with shape: broadcast(img1.shape[:-3], img2.shape[:-3]).
  """
   _, _, checks = _verify_compatible_image_shapes(img1, img2)
  with ops.control_dependencies(checks):
    img1 = array_ops.identity(img1)

  # Need to convert the images to float32.  Scale max_val accordingly so that
  # SSIM is computed correctly.
  max_val = math_ops.cast(max_val, img1.dtype)
  max_val = convert_image_dtype(max_val, dtypes.float32)
  img1 = convert_image_dtype(img1, dtypes.float32)
  img2 = convert_image_dtype(img2, dtypes.float32)
  ssim_per_channel, _ = _ssim_per_channel(img1, img2, max_val)
  # Compute average over color channels.
  return math_ops.reduce_mean(ssim_per_channel, [-1])

2、Python实现

该实现摘抄自：https://blog.csdn.net/sinat_36438332/article/details/88173349

作者实际上采用σ=1.5的高斯加权函数计算每个图像块(图像块大小11x11)的均值和标准差，目的是防止undesirable “blocking” artifacts：

生成高斯核的函数参考了RaymondMG的博客
由于使用的高斯函数圆对称，因此相关操作和卷积操作结果相同，这里只定义了相关操作，卷积在相关操作前将核旋转180°即可。

import cv2
import numpy as np
import time
from numba import jit,njit

#相关操作
#由于使用的高斯函数圆对称，因此相关操作和卷积操作结果相同
@njit
def correlation(img,kernal):
    kernal_heigh = kernal.shape[0]
    kernal_width = kernal.shape[1]
    cor_heigh = img.shape[0] - kernal_heigh + 1
    cor_width = img.shape[1] - kernal_width + 1
    result = np.zeros((cor_heigh, cor_width), dtype=np.float64)
    for i in range(cor_heigh):
        for j in range(cor_width):
            result[i][j] = (img[i:i + kernal_heigh, j:j + kernal_width] * kernal).sum()
    return result

#产生二维高斯核函数
#这个函数参考自：https://blog.csdn.net/qq_16013649/article/details/78784791
@jit
def gaussian_2d_kernel(kernel_size=11, sigma=1.5):
    kernel = np.zeros([kernel_size, kernel_size])
    center = kernel_size // 2

    if sigma == 0:
        sigma = ((kernel_size - 1) * 0.5 - 1) * 0.3 + 0.8

    s = 2 * (sigma ** 2)
    sum_val = 0
    for i in range(0, kernel_size):
        for j in range(0, kernel_size):
            x = i - center
            y = j - center
            kernel[i, j] = np.exp(-(x ** 2 + y ** 2) / s)
            sum_val += kernel[i, j]
    sum_val = 1 / sum_val
    return kernel * sum_val


#ssim模型
@jit
def ssim(distorted_image,original_image,window_size=11,gaussian_sigma=1.5,K1=0.01,K2=0.03,alfa=1,beta=1,gama=1):
    distorted_image=np.array(distorted_image,dtype=np.float64)
    original_image=np.array(original_image,dtype=np.float64)
    if not distorted_image.shape == original_image.shape:
        raise ValueError("Input Imagees must has the same size")
    if len(distorted_image.shape) > 2:
        raise ValueError("Please input the images with 1 channel")
    kernal=gaussian_2d_kernel(window_size,gaussian_sigma)

    #求ux uy ux*uy ux^2 uy^2 sigma_x^2 sigma_y^2 sigma_xy等中间变量
    ux=correlation(distorted_image,kernal)
    uy=correlation(original_image,kernal)
    distorted_image_sqr=distorted_image**2
    original_image_sqr=original_image**2
    dis_mult_ori=distorted_image*original_image
    uxx=correlation(distorted_image_sqr,kernal)
    uyy=correlation(original_image_sqr,kernal)
    uxy=correlation(dis_mult_ori,kernal)
    ux_sqr=ux**2
    uy_sqr=uy**2
    uxuy=ux*uy
    sx_sqr=uxx-ux_sqr
    sy_sqr=uyy-uy_sqr
    sxy=uxy-uxuy
    C1=(K1*255)**2
    C2=(K2*255)**2
    #常用情况的SSIM
    if(alfa==1 and beta==1 and gama==1):
        ssim=(2*uxuy+C1)*(2*sxy+C2)/(ux_sqr+uy_sqr+C1)/(sx_sqr+sy_sqr+C2)
        return np.mean(ssim)
    #计算亮度相似性
    l=(2*uxuy+C1)/(ux_sqr+uy_sqr+C1)
    l=l**alfa
    #计算对比度相似性
    sxsy=np.sqrt(sx_sqr)*np.sqrt(sy_sqr)
    c=(2*sxsy+C2)/(sx_sqr+sy_sqr+C2)
    c=c**beta
    #计算结构相似性
    C3=0.5*C2
    s=(sxy+C3)/(sxsy+C3)
    s=s**gama
    ssim=l*c*s
    return np.mean(ssim)

另外，代码中的@njit和@jit可以注释掉，我在这里是为了加速代码运行(实际加速效果不大)。

参考博客：

1、ssim算法原理

2、SSIM(结构相似性)-数学公式及python实现

3、图像质量评估算法 SSIM(结构相似性)

4、结构相似形特征(SSIM)原理简介及python实现

5、https://baike.baidu.com/item/SSIM/2091025?fr=aladdin