图像无参考评价模型BRISQUE原理简介以及Python(Pytorch)实现

本文主要是介绍图像无参考评价模型BRISQUE原理简介以及Python(Pytorch)实现，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

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{aligned} &\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{aligned}$
$C = 1$ ，是为了防止分母为0。
可见这其实很类似数据预处理过程中的z-score归一化过程。
归一化之后的图像数据不再用 $I$ 表示,而是用 $M S C N$ (mean subtracted contrast normalized coefficients)表示。
局部亮度归一化的一种生物学解释是，人来在感知失真时，会存在一种“对比掩蔽效应”。这种效应通俗的解释就是，在图像纹理很丰富的地方，图像的失真往往很难被感受到，复杂的纹理掩盖了失真。这种归一化的过程，在其他的论文里面也称为Divisive Normalization。

2.利用广义高斯分布拟合 $M S C N$ 获得特征

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

零均值广义高斯分布(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} d t \quad a>0$
利用GGD拟合分布获得的形状参数 $α, σ^2)$ 作为特征f1和特征f2。关于广义高斯分布的参数拟合方法实现问题，我在另一篇博客单独描述，这里主要介绍IQA模型原理，不再叙述。

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

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

$f\left(x ; v, \sigma_{l}^{2}, \sigma_{r}^{2}\right)=\left\{\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 \geq 0 \end{array}\right.$
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拟合获得的形状参数 $η, ν, σ_l^2 , σ_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 // 2if sigma == 0:sigma = ((kernel_size - 1) * 0.5 - 1) * 0.3 + 0.8s = 2 * (sigma ** 2)sum_val = 0for i in range(0, kernel_size):for j in range(0, kernel_size):x = i - centery = j - centerkernel[i, j] = math.exp(-(x ** 2 + y ** 2) / s)sum_val += kernel[i, j]sum_val = 1 / sum_valreturn kernel * sum_valdef 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_sqdef 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_ru2=np.mean(vec**2)m1=np.mean(np.abs(vec))r_=m1**2/u2R_=r_*(gamma_**3+1)*(gamma_+1)/((gamma_**2+1)**2)diff=(R_-r_alpha)**2alpha_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*const2return alpha_param,eta,sigma_l**2,sigma_r**2def brisque_feature(imdist:Type[Union[torch.Tensor,np.ndarray]],device='cuda'):#算法需要输入为灰度图像，像素值0-255if type(imdist)==np.ndarray:assert imdist.ndim==2 or imdist.ndim==3if 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()==4assert imdist.shape[1]==1 or imdist.shape[1]==3if torch.max(imdist)<=1:imdist = imdist * 255# RGB to Grayif imdist.shape[1]==3:imdist=imdist[:,0,:]*0.299+imdist[:,1,:]*0.587+imdist[:,2,:]*0.114# GPU is much much fasterif 'cuda' in device:imdist=imdist.half().to(device)elif device=='cpu':imdist=imdist.float().to(device)else:raise ValueError('cpu or cuda',device)# 算法主体scalenum = 2window=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*musigma=torch.sqrt(torch.abs(F.conv2d(imdist*imdist,window,stride=1,padding=3)-mu_sq))structdis = (imdist - mu) / (sigma + 1)del mu, mu_sq,sigmafeat[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_structdispair=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)