结构相似性(SSIM)原理及其实现

一、结构相似性

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

二、SSIM指数

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

SSIM测量系统由三个对比模块组成：亮度、对比度、结构

2.1 亮度对比函数

将图像的平均灰度作为亮度测量的估计
$${\mu }_X=\frac{1}{H\times M}\sum _{i=1}^H\sum _{j=1}^{M}X(i,j)$$
于是两幅图像的亮度对比函数
$$l(x,y)=\frac{2{\mu }_x{\mu }_y+C_1}{{\mu }_x^2+{\mu }_y^2+C_1}$$

2.2 对比度对比函数

将图像的标准差作为对比度测量估计
$${\sigma }_X=(\frac{1}{H+W-1}\sum _{i=1}^H\sum _{j=1}^M(X(i,j)-{\mu }_X{)}^2{)}^{\frac{1}{2}}$$
于是两幅图像的对比度对比函数
$$c(x,y)=\frac{2{\sigma }_x{\sigma }_y+C_2}{{\sigma }_x^2+{\sigma }_y^2+C_2}$$

2.3 结构对比函数

结构对比函数
$$s(x,y)=\frac{{\sigma }_{xy}+C_3}{{\sigma }_x{\sigma }_y+C_3}$$

2.4 SSIM测量函数

最后， 以上三个对比函数构成了SSIM函数
$$\begin{array}{rl}SSIM(x,y)& =f(l(x,y),c(x,y),s(x,y))\\ & =[l(x,y){]}^{\alpha }[c(x,y){]}^{\beta }[s(x,y){]}^{\gamma }\end{array}$$
其中， $\alpha 、\beta 、\gamma >0$ ，用来调整这三个模块的重要性
假设$\alpha 、\beta 、\gamma 都为1$， $C_3=C_2/2$ ,则
$$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函数的值域为[0, 1], 值越大说明图像失真越小，两幅图像越相似

2.5 SSIM函数满足的三个条件

• 对称性：$S(x,y)=S(y,x)$
• 有界性：$S(x,y)\le 1$
• 最大值唯一性：当且仅当x=y时，$S(x,y)=1$

三、MSSIM

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

可以利用滑动窗将图像分块，令分块总数为N，考虑到窗口形状对分块的影响，采用高斯加权计算每一窗口的均值、方差以及协方差，然后计算对应块的结构相似度SSIM，最后将平均值作为两图像的结构相似性度量，即平均结构相似性MSSIM：

四、实现

import torch
import torch.nn.functional as F
from math import exp
import numpy as np

# 创建一维高斯分布向量
def gaussian(window_size, sigma):
    gauss = torch.Tensor([exp(-(x-window_size//2)**2)/float(2*sigma**2) for x in range(window_size)])
    return gauss/gauss.sum()

# 创建高斯核
def create_window(window_size, channel=1):
    # unqueeze(1) 在第二维上增加一个维度
    _1D_window = gaussian(window_size, 1.5).unsqueeze(1)
    # t() 转置
    _2D_window = _1D_window.mm(_1D_window.t()).float().unsqueeze(0).unsqueeze(0)
    window = _2D_window.expand(channel, 1, window_size, window_size).contiguous()
    return window


# 计算ssim
# 采用归一化的高斯核来 代替计算像素的均值
def ssim(img1, img2, window, window_size, channel=1, size_average=True):
    
    (
    mu1 = F.conv2d(img1, window, padding=window_size//2, stride=1, groups=channel)
    mu2 = F.conv2d(img2, window, padding=window_size//2, stride=1, groups=channel)
    
    mu1_sq = mu1.pow(2)
    mu2_sq = mu2.pow(2)
    mu1_mu2 = mu1 *mu2
    
    sigma1_sq = F.conv2d(img1*img1, window, padding=window_size//2, stride=1, groups=channel) - mu1_sq
    sigma2_sq = F.conv2d(img2*img2, window, padding=window_size//2, stride=1, groups=channel) - mu2_sq
    sigma_12 = F.conv2d(img1*img2, window, padding=window_size//2, stride=1, groups=channel) - mu1_mu2
    
    c1 = 0.01 **2
    c2 = 0.03 **2
    
    ssim_map = (2*(mu1_mu2 +c1)*(2*sigma_12 + c2)) / ((mu1_sq + mu2 + c1)*(sigma1_sq + sigma2_sq + c2))
    
    if size_average:
        return ssim_map.mean()
    else:
        return ssim_map.mean(1).mean(1).mean(1)
    
    
class SSIM(torch.nn.Module):
    def __init__(self, window_size=11, channel, size_average=True):
        super(SSIM, self).__init__()
        self.window_size = window_size
        self.size_average = size_average
        self.channel = channel
        self.window = create_window(window_size, channel)
        
    def forward(self, img1, img2):
        (_, channel, _, _) = img1.size()
        
        if channel == self.channel and self.window.data.type() == img1.data.type():
            window = self.window
        else:
            window = create_window(self.window_size, channel)
            if img1.is_cuda:
                window.cuda(img1.get_device())
            window = window.type_as(img1)
            self.window = window
            self.channel = channel
        
        return ssim(img1, img2, self.window, self.window_size, channel, self.size_average)