图像全参考评价之梯度幅度相似性偏差(GMSD)原理简介以及Python(Pytorch)实现

梯度幅度相似性偏差(Gradient Magnitude Similarity Deviation)是2014年zhang lei等人在论文《Gradient Magnitude Similarity Deviation: A Highly Efficient Perceptual Image Quality Index》提出的一种图像全参考评价(FR-IQA)方法，具有准确度高、计算量少的特点。

原理简介

梯度幅值能够反映结构信息，仅使用梯度幅值作为特征就能产生准确度较高的图片质量预测分数。自然图像往往有着丰富类型的局部结构，不同的结构在失真中会有不同的梯度幅值退化。基于局部质量退化在图像全局上的变化可以反映图像的质量这一观点，文章提出了计算局部梯度幅值相似性来衡量局部图像质量，计算局部图像质量的标准差来衡量图像全局的质量。

1.计算局部梯度幅值相似性

作者使用Prewitt算子计算图像梯度，使用其他算子如 Sobel and Scharr 得到的效果相似。

${\mathbf{h}}_x=\left[\begin{array}{lll}1/3& 0& -1/3\\ 1/3& 0& -1/3\\ 1/3& 0& -1/3\end{array}\right],{\mathbf{h}}_y=\left[\begin{array}{ccc}1/3& 1/3& 1/3\\ 0& 0& 0\\ -1/3& -1/3& -1/3\end{array}\right]$

然后计算梯度幅值 ${\mathbf{m}}_r$、 ${\mathbf{m}}_d$：

${\mathbf{m}}_r(i)=\sqrt{{\left(\mathbf{r}\otimes {\mathbf{h}}_x\right)}^2(i)+{\left(\mathbf{r}\otimes {\mathbf{h}}_y\right)}^2(i)}$
${\mathbf{m}}_d(i)=\sqrt{{\left(\mathbf{d}\otimes {\mathbf{h}}_x\right)}^2(i)+{\left(\mathbf{d}\otimes {\mathbf{h}}_y\right)}^2(i)}$

梯度幅值相似性(GMS)如下：

$GMS(i)=\frac{2{\mathbf{m}}_r(i){\mathbf{m}}_d(i)+c}{{\mathbf{m}}_r^2(i)+{\mathbf{m}}_d^2(i)+c}$

c是一个常数，作用是为了防止分母为0。

2.标准差池化

SSIM等方法使用的是平均池化方法，这种方法计算简便，认为图像每块区域的质量对全局质量的贡献相同。也有研究人员提出加权池化的方法，但是这种方法对效果提升作用不大且提高了计算量。自然图像具有丰富的结构，不同结构对于失真的响应不同，如模糊在“平坦”部位的影响就比在“纹理”部位的影响小。基于局部质量退化在图像全局上的变化可以反映图像的质量这一理念，文章提出了 Standard Deviation Pooling。

$GMSD=\sqrt{\frac{1}{N}{\sum }_{i=1}^N(GMS(i)-GMSM{)}^2}$

代码实现

在实验应用中，文章首先使用2x2均值滤波器分别对原始图像和失真图像进行滤波，然后进行间隔为2的下采样。

论文的MATLAB源码如下：

function [score, quality_map] = GMSD(Y1, Y2)

% GMSD - measure the image quality of distorted image 'Y2' with the reference image 'Y1'.
% 
% inputs:
% 
% Y1 - the reference image (grayscale image, double type, 0~255)
% Y2 - the distorted image (grayscale image, double type, 0~255)
% 
% outputs:

% score: distortion degree of the distorted image
% quality_map: local quality map of the distorted image

% This is an implementation of the following algorithm:
% Wufeng Xue, Lei Zhang, Xuanqin Mou, and Alan C. Bovik, 
% "Gradient Magnitude Similarity Deviation: A Highly Efficient Perceptual Image Quality Index",
% http://www.comp.polyu.edu.hk/~cslzhang/IQA/GMSD/GMSD.htm

T = 170; 
Down_step = 2;
dx = [1 0 -1; 1 0 -1; 1 0 -1]/3;
dy = dx';

aveKernel = fspecial('average',2);
aveY1 = conv2(Y1, aveKernel,'same');
aveY2 = conv2(Y2, aveKernel,'same');
Y1 = aveY1(1:Down_step:end,1:Down_step:end);
Y2 = aveY2(1:Down_step:end,1:Down_step:end);

IxY1 = conv2(Y1, dx, 'same');     
IyY1 = conv2(Y1, dy, 'same');    
gradientMap1 = sqrt(IxY1.^2 + IyY1.^2);

IxY2 = conv2(Y2, dx, 'same');     
IyY2 = conv2(Y2, dy, 'same');

gradientMap2 = sqrt(IxY2.^2 + IyY2.^2);
quality_map = (2*gradientMap1.*gradientMap2 + T) ./(gradientMap1.^2+gradientMap2.^2 + T);
score = std2(quality_map);

python实现如下

由于涉及到卷积操作，我基于Pytorch来实现，Pytorch计算卷积比较快，并且可以通过GPU进行加速。

import torch
import torch.nn.functional as F
import skimage.io
import numpy as np
from typing import Type, Any, Callable, Union, List, Optional

def gmsd(dis_img:Type[Union[torch.Tensor,np.ndarray]],ref_img:Type[Union[torch.Tensor,np.ndarray]],c=170,device='cuda'):
    # 输入类型检查
    if type(dis_img) == np.ndarray:
        assert dis_img.ndim == 2 or dis_img.ndim == 3
        if dis_img.ndim == 2:
            dis_img = torch.from_numpy(dis_img).unsqueeze(0).unsqueeze(0)
        else:
            dis_img = torch.from_numpy(dis_img).unsqueeze(0)

    if type(ref_img) == np.ndarray:
        assert ref_img.ndim == 2 or ref_img.ndim == 3
        if ref_img.ndim == 2:
            ref_img = torch.from_numpy(ref_img).unsqueeze(0).unsqueeze(0)
        else:
            ref_img = torch.from_numpy(ref_img).unsqueeze(0)
    # 算法需要输入为灰度图像，像素值0-255
    if torch.max(dis_img) <= 1:
        dis_img = dis_img * 255
    if torch.max(ref_img) <= 1:
        ref_img = ref_img * 255

    '''算法主体'''
    hx=torch.tensor([[1/3,0,-1/3]]*3,dtype=torch.float).unsqueeze(0).unsqueeze(0).to(device)#Prewitt算子
    ave_filter=torch.tensor([[0.25,0.25],[0.25,0.25]],dtype=torch.float).unsqueeze(0).unsqueeze(0).to(device)#均值滤波核
    down_step=2#下采样间隔
    hy=hx.transpose(2,3)

    dis_img=dis_img.float().to(device)
    ref_img=ref_img.float().to(device)

    #均值滤波
    ave_dis=F.conv2d(dis_img,ave_filter,stride=1)
    ave_ref=F.conv2d(ref_img,ave_filter,stride=1)
    #下采样
    ave_dis_down=ave_dis[:,:,0::down_step,0::down_step]
    ave_ref_down=ave_ref[:,:,0::down_step,0::down_step]
    #计算mr md等中间变量
    mr_sq=F.conv2d(ave_ref_down,hx)**2+F.conv2d(ave_ref_down,hy)**2
    md_sq=F.conv2d(ave_dis_down,hx)**2+F.conv2d(ave_dis_down,hy)**2
    mr=torch.sqrt(mr_sq)
    md=torch.sqrt(md_sq)
    GMS=(2*mr*md+c)/(mr_sq+md_sq+c)
    GMSD=torch.std(GMS.view(-1))
    return GMSD

if __name__=='__main__':
    img_ref = skimage.io.imread(r'D:\DataBase\IMAGE\databaserelease2\refimgs\bikes.bmp', as_gray=True)
    img_dis=skimage.io.imread(r'D:\DataBase\IMAGE\databaserelease2\jp2k\img225.bmp', as_gray=True)
    score=gmsd(img_dis,img_ref)
    print(score)