0 引言

 自然场景统计(NSS)经常被应用于图像的质量评价。人们发现，用高质量设备采集的自然图像(自然场景)有着一定的统计特征(如服从一些类高斯分布)，而图像的失真会使这些统计特征发生改变，因此利用图像的一些统计特征作为图像特征，可以完成图像的无参考评价。
BRISQUE来自于论文《No-Reference Image Quality Assessment in the Spatial Domain》，是一个经典的利用NSS进行NR-IQA的模型。观察到自然图像局部归一化亮度系数(MSCN)，强烈趋向于单位正态高斯特性，因此作者假设失真会改变MSCN的分布状况，并基于此提取特征。
 基于自然场景统计的工作还有很多，如BLIINDS-II、VIF、DIIVINE等，但是这些工作往往需要对图形进行分解(DCT、小波等)。BRISQUE的优美之处在于，它不需要对图像进行频域分解，仅仅使用简单的归一化过程，就使得数据呈现有规律的分布。模型简单且高效，可扩展性强，计算的复杂度较低。
 后来经过实验发现，四个方向内积的AGGD特征是模型的主要特征，他们给模型的性能带来主要的贡献。这表明，图像的结构是影响人类对图像质量感知的主要因素。在质量评价中，绝大部分模型都会考虑到图像的结构信息，如经典的SSIM，或者LGP、GMSD、GM-LOG使用梯度信息模型，以及其他一些频域分解的模型。而BRISQUE使用四个方向的点积来表示图像的结构，这种想法十分巧妙。

1 原理简介

1.局部亮度归一化

 如上图左列所示，像素与其周围的像素具有很高的相关性，因为图像函数通常是分段平滑的。我们在左列可以观察到一种对角结构。图像的归一化能够大大减少像素的相关性，如右列所示。
局部亮度归一化过程如下式：
$$\hat{I}(i,j)=\frac{I(i,j)-\mu (i,j)}{\sigma (i,j)+C}$$
其中：
$$\begin{array}{rl}& \mu (i,j)=\sum _{k=-K}^K\sum _{l=-L}^L{w}_{k,l}{I}_{k,l}(i,j)\\ & \sigma (i,j)=\sqrt{\sum _{k=-K}^K\sum _{l=-L}^L{w}_{k,l}{\left(I_{k,l}(i,j)-\mu (i,j)\right)}^2}\end{array}$$
 $C=1$，是为了防止分母为0。
 可见这其实很类似数据预处理过程中的z-score归一化过程。
 归一化之后的图像数据不再用$I$表示,而是用$MSCN$(mean subtracted contrast normalized coefficients)表示。
 局部亮度归一化的一种生物学解释是，人来在感知失真时，会存在一种“对比掩蔽效应”。这种效应通俗的解释就是，在图像纹理很丰富的地方，图像的失真往往很难被感受到，复杂的纹理掩盖了失真。这种归一化的过程，在其他的论文里面也称为Divisive Normalization。

2.利用广义高斯分布拟合$MSCN$获得特征

  自然图像(未失真)的$MSCN$具有一定的统计特征(如类高斯分布)，而失真改变这一分布。下图可以观察到不同失真类型图像和自然图像具有不同分布特征，文献使用广义高斯分布(GGD)对图像$MSCN$的分布进行拟合。

零均值广义高斯分布(GGD)的定义如下：
$$f\left(x;\alpha ,{\sigma }^2\right)=\frac{\alpha }{2\beta \Gamma (1/\alpha )}\exp \left(-{\left(\frac{|x|}{\beta }\right)}^{\alpha }\right)$$
where
$$\beta =\sigma \sqrt{\frac{\Gamma (1/\alpha )}{\Gamma (3/\alpha )}}$$
and $\Gamma (\cdot )$ is the gamma function:
$$\Gamma (a)={\int }_0^{\infty }{t}^{a-1}{e}^{-t}dta>0$$
 利用GGD拟合分布获得的形状参数$(\alpha ,{\sigma }^2)$作为特征f1和特征f2。关于广义高斯分布的参数拟合方法实现问题，我在另一篇博客单独描述，这里主要介绍IQA模型原理，不再叙述。

3.利用非对称广义高斯分布拟合MSCN相邻系数内积

 如前面所说，图像的结构信息具有很大的价值，因此论文也研究了相邻MSCN系数之间的统计特征来自作为图像的结构信息。论文在四个不同方向计算相邻系数内积，并利用非对称广义高斯分布(AGGD)拟合分布特征。四个方向的内积定义如下：
$$\begin{array}{l}H(i,j)=\hat{I}(i,j)\ast \hat{I}(i,j+1)\\ V(i,j)=\hat{I}(i,j)\ast \hat{I}(i+1,j)\\ D1(i,j)=\hat{I}(i,j)\ast \hat{I}(i+1,j+1)\\ D2(i,j)=\hat{I}(i,j)\ast \hat{I}(i+1,j-1)\end{array}$$
 下图展示了其中一个方向上内积数据的分布情况，可以发现失真的引入会使得分布带来显著改变。这个分布，论文使用非对称广义高斯分布AGGD拟合。

零均值非对称广义高斯分布(AGGD)的定义如下：

$$f\left(x;v,{\sigma }_l^2,{\sigma }_r^2\right)=\{\begin{array}{ll}\frac{v}{\left({\beta }_1+{\beta }_r\right)\Gamma \left(\frac{1}{v}\right)}\exp \left(-{\left(\frac{-x}{{\beta }_l}\right)}^v\right)& x<0\\ \frac{v}{\left({\beta }_l+{\beta }_r\right)\Gamma \left(\frac{1}{v}\right)}\exp \left(-{\left(\frac{x}{{\beta }_r}\right)}^v\right)& x\ge 0\end{array}$$
where
$$\begin{array}{l}{\beta }_l={\sigma }_l\sqrt{\frac{\Gamma \left(\frac{1}{v}\right)}{\Gamma \left(\frac{3}{v}\right)}}\\ {\beta }_r={\sigma }_r\sqrt{\frac{\Gamma \left(\frac{1}{v}\right)}{\Gamma \left(\frac{3}{v}\right)}}\end{array}$$
$$\eta =\left({\beta }_r-{\beta }_l\right)\frac{\Gamma \left(\frac{2}{v}\right)}{\Gamma \left(\frac{1}{v}\right)}$$

 利用AGGD拟合获得的形状参数$(\eta ,\nu ,{\sigma }_l^2,{\sigma }_r^2)$将作为特征f3-f18(4个方向 x 4个特征)。关于非对称广义高斯分布的参数拟合方法实现问题，我亦在另一篇博客单独描述，这里主要介绍IQA模型原理，不再叙述。

 至此，对图像的18个特征提取就完成了，具体的特征内容如下表。

 需要指出的是，多尺度的特征经常在IQA模型中使用。论文采用了多尺度特征的形式，具体操作为首先对原始图像采集f1-f18这18维特征，然后对图像进行一次下采样，在下采样的图像上再提取一次这18维的特征，这样一共获得36维特征。

2 代码实现

• 文献的作者是给出MATLAB代码的，地址如下：
  http://live.ece.utexas.edu/research/quality/BRISQUE_release.zip
  一般用迅雷就能下下来。我也存到网盘上了，下不了的可以用网盘：
  链接：https://pan.baidu.com/s/1RxlaOAlGGMLWZV7X7sTqOQ
  提取码：bab3
• 最近发现作者也更新了一个C++版本的代码：
  http://live.ece.utexas.edu/research/Quality/brisque_revised.zip
  网盘链接：https://pan.baidu.com/s/1Ei204UN3xBWMTWS0Zakkig
  提取码：zt71
• 本文的python实现：https://github.com/guyueyao/torchBRISUQE

本文主要给出BRISQUE算法的特征提取部分：

import torch
import torch.nn.functional as F
import math
from scipy.special import gamma
import numpy as np
from typing import Type, Any, Callable, Union, List, Optional
def gaussian_2d_kernel(kernel_size, sigma):
    kernel = torch.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] = math.exp(-(x ** 2 + y ** 2) / s)
            sum_val += kernel[i, j]
    sum_val = 1 / sum_val
    return kernel * sum_val

def estimate_GGD_parameters(vec):
    vec=vec.view(-1)
    gam =np.arange(0.2,10.0,0.001)
    r_gam = (gamma(1/gam)*gamma(3/gam))/((gamma(2/gam))**2)#根据候选的γ计算r(γ)
    sigma_sq=torch.mean((vec)**2).to('cpu').numpy()
    E=torch.mean(torch.abs(vec)).to('cpu').numpy()
    r=sigma_sq/(E**2)#根据sigma^2和E计算r(γ)
    diff=np.abs(r-r_gam)
    gamma_param=gam[np.argmin(diff, axis=0)]
    return gamma_param,sigma_sq

def estimate_AGGD_parameters(vec):
    vec = vec.to('cpu').numpy()
    alpha =np.arange(0.2,10.0,0.001)#产生候选的α
    r_alpha=((gamma(2/alpha))**2)/(gamma(1/alpha)*gamma(3/alpha))#根据候选的γ计算r(α)
    sigma_l=np.sqrt(np.mean(vec[vec<0]**2))
    sigma_r=np.sqrt(np.mean(vec[vec>0]**2))
    gamma_=sigma_l/sigma_r
    u2=np.mean(vec**2)
    m1=np.mean(np.abs(vec))
    r_=m1**2/u2
    R_=r_*(gamma_**3+1)*(gamma_+1)/((gamma_**2+1)**2)
    diff=(R_-r_alpha)**2
    alpha_param=alpha[np.argmin(diff, axis=0)]
    const1 = np.sqrt(gamma(1 / alpha_param) / gamma(3 / alpha_param))
    const2 = gamma(2 / alpha_param) / gamma(1 / alpha_param)
    eta =(sigma_r-sigma_l)*const1*const2
    return alpha_param,eta,sigma_l**2,sigma_r**2


def brisque_feature(imdist:Type[Union[torch.Tensor,np.ndarray]],device='cuda'):
    
    #算法需要输入为灰度图像，像素值0-255
    if type(imdist)==np.ndarray:
        assert imdist.ndim==2 or imdist.ndim==3
        if imdist.ndim==2:
            imdist=torch.from_numpy(imdist).unsqueeze(0).unsqueeze(0)
        else:
            imdist = torch.from_numpy(imdist).unsqueeze(0)
    # input (Batch,1,H,W)
    assert imdist.dim()==4
    assert imdist.shape[1]==1 or imdist.shape[1]==3
    
    if torch.max(imdist)<=1:
        imdist = imdist * 255
    # RGB to Gray
    if imdist.shape[1]==3:
        imdist=imdist[:,0,:]*0.299+imdist[:,1,:]*0.587+imdist[:,2,:]*0.114
    # GPU is much much faster
    if 'cuda' in device:
        imdist=imdist.half().to(device)
    elif device=='cpu':
        imdist=imdist.float().to(device)
    else:
        raise ValueError('cpu or cuda',device)

    # 算法主体
    scalenum = 2
    window=gaussian_2d_kernel(7,7/6).unsqueeze(0).unsqueeze(0).float().to(device)
    if 'cuda' in device:
        window=window.half()

    feat=np.zeros((18*scalenum,))
    for i in range(scalenum):
        mu=F.conv2d(imdist,window,stride=1,padding=3)
        mu_sq=mu*mu
        sigma=torch.sqrt(torch.abs(F.conv2d(imdist*imdist,window,stride=1,padding=3)-mu_sq))
        structdis = (imdist - mu) / (sigma + 1)
        del mu, mu_sq,sigma
        feat[i*18],feat[i*18+1] = estimate_GGD_parameters(structdis)

        shifts = [(0,1),(1,0),(1,1),(-1,1)]
        for itr_shift in range(4):
            shifted_structdis=structdis.roll(shifts[itr_shift],(2,3))
            pair=structdis*shifted_structdis
            pair=pair.view(-1)
            feat[i*18+2+itr_shift*4],feat[i*18+3+itr_shift*4],feat[i*18+4+itr_shift*4],feat[i*18+5+itr_shift*4]=estimate_AGGD_parameters(pair)

        imdist=F.interpolate(imdist,scale_factor=(0.5,0.5),mode='bilinear')
    return feat

算法主要是利用Pytorch实现，在GPU上跑的话相对较快。上述代码仅为特征提取部分，模型训练时需要提取训练集图像特征，之后利用特征和对应的质量分数训练SVR；测试时需要先提取测试图像特征，之后应用SVR回归。下面的不完整代码简单描述了这个过程：

from sklearn.svm import SVR
from brisque_feature import brisque_feature
import numpy as np
import skimage.io
# 训练
X=[]
Y=[]
for img_name,img_dmos in train_set:
	img=skimage.io.imread(os.path.join(root_path,img_name),as_gray=True)#读图像
    feat=brisque_feature(img)#提特征
    X.append(feat)
    Y.append(img_dmos)
X=np.array(X)
Y=np.array(Y)
svr = SVR(kernel='rbf', C=C, gamma=Gamma)#训练SVR
svr.fit(X,Y)

# 测试
img=skimage.io.imread(os.path.join(test_img_path),as_gray=True)
feat=brisque_feature(img)
predicted_score=svr.predict(feat)

默认使用GPU进行计算，如果没有GPU可以选择使用CPU：

feat=brisque_feature(img,'cpu')