NIQE 概述
NIQE（Natural Image Quality Evaluator，自然图像质量评估器）是一种无参考图像质量评价算法，用于评估图像的自然度和失真程度。该算法旨在解决传统的基于参考图像的评估指标（如峰值信噪比PSNR、结构相似性SSIM等）在评价图像质量时与人类视觉感知不一致的问题。
NIQE算法的核心思想是利用自然图像的统计特征来建立一个模型，用于估计待评价图像的质量。具体实现过程包括以下几个步骤：
- 图像预处理：将待评价图像转换为灰度图像，并调整图像大小和对比度等。
- 特征提取：从预处理后的图像中提取一系列特征，包括局部统计特征和全局统计特征。局部统计特征通常包括图像的梯度、对比度和结构信息等，全局统计特征则包括图像的均值、方差和灰度共生矩阵等。
- 模型建立：基于提取的特征，使用机器学习方法（如支持向量机SVM或随机森林Random Forest）来构建一个多元高斯模型。
- 质量评估：将给定的测试图像的NIQE指标表示为从测试图像中提取的NSS特征的MVG模型与从自然图像语料中提取的质量感知特征的MVG模型之间的距离。这个距离反映了测试图像与自然图像的统计差异，从而可以用来评估图像的质量。
NIQE算法的优点在于它不需要参考原始的高质量图像，因此在实际应用中更加灵活和方便。此外，由于其与人类视觉感知的一致性较高，NIQE算法在图像质量评估领域得到了广泛应用。
需要注意的是，虽然NIQE算法在图像质量评估方面具有较高的准确性和可靠性，但它仍然受到一些限制。例如，对于某些特定类型的失真（如噪声、模糊等），NIQE算法可能无法给出准确的评估结果。因此，在实际应用中，需要根据具体情况选择合适的图像质量评估算法。
一、NIQE 测试
1. 生成多等级扭曲图像代码
ImageList = dir("../Capture/*.jpg");
PicNum = size(ImageList);
ImagePath = strcat(strcat(ImageList(PicNum(1) - 1).folder, '\'), ImageList(PicNum(1) - 1).name);
I = imread(ImagePath);
I_gray = rgb2gray(I);
figure,
imshow(I_gray)
figure,
% I_Trans = imnoise(I_gray,"gaussian", 0, 0.21); % Guassian Transform for Origin Image File. [0.01, 0.05, 0.09, 0.13, 0.17
% I_Trans = imnoise(I_gray,'salt & pepper', 0.17);
% I_Trans = imnoise(I_gray,'speckle',0.01);
% I_Trans = imadjust(I_gray,[0.1 0.9],[0.7 0.9]);
% h = fspecial('motion', 50, 50); I_Trans = imfilter(I_gray, h);
h = fspecial('gaussian', [7 7], 2); I_Trans = imfilter(I_gray, h, 'replicate');
imshow(I_Trans)
ImageNoiseTransformLevel = [0.01, 0.05, 0.09, 0.13, 0.17];
ImageDarkTransformLevel = [0.5, 0.4, 0.3, 0.2, 0.1];
ImageBrightenTransformLevel = [0.4, 0.5, 0.6, 0.7, 0.8];
ImageMotionBlurTransformLevel = [10, 20, 30, 40, 50];
ImageGaussianBlurTransformLevel = [1, 5, 10, 15, 20];
for Lev = 1:5
for i = 1:PicNum(1) % 3
ImagePath = strcat(strcat(ImageList(i).folder, '\'), ImageList(i).name);
I = imread(ImagePath);
try
I_gray = rgb2gray(I);
catch ME
I_gray = I;
end
% I_Trans = imnoise(I_gray,"gaussian", 0, ImageNoiseTransformLevel); % [0.01, 0.05, 0.09, 0.13, 0.17]
% I_Trans = imnoise(I_gray,'speckle', ImageNoiseTransformLevel(Lev)); % [0.01, 0.05, 0.09, 0.13, 0.17]
% I_Trans = imnoise(I_gray,'salt & pepper', ImageNoiseTransformLevel(Lev));
% I_Trans = imadjust(I_gray,[0.1 0.9],[0.1 ImageDarkTransformLevel(Lev)]);
% I_Trans = imadjust(I_gray,[0.1 0.9],[ImageBrightenTransformLevel(Lev) 0.9]);
% h = fspecial('motion', ImageMotionBlurTransformLevel(Lev), ImageMotionBlurTransformLevel(Lev)); I_Trans = imfilter(I_gray, h);
h = fspecial('gaussian', [7 7], ImageGaussianBlurTransformLevel(Lev)); I_Trans = imfilter(I_gray, h, 'replicate');
imshow(I_Trans)
SavePath = strcat(strcat(strcat(ImageList(i).folder, '\GaussianBlur\L'), num2str(Lev)), ImageList(i).name);
imwrite(I_Trans, SavePath);
end
end
🌻 生成图像结果：
2. 使用NIQE进行拟合得分
% 使用自定义特征模型计算NIQE分数, 将一组自然图像加载到一个图像数据存储中
setDir = ['Z:\LT100\MaMiao\AWorkSpace\项目计划\LT-200\lt-200develop' ...
'\01_可行性研究\OpticalDesign\ImageQuanlity\ImageResource\Capture\FittingImage'];
imds = imageDatastore(setDir,'FileExtensions',{'.jpg'});
% 使用图像数据存储训练自定义NIQE模型
model = fitniqe(imds);
img_folder = ['Z:\LT100\MaMiao\AWorkSpace\项目计划\LT-200\lt-200develop' ...
'\01_可行性研究\OpticalDesign\ImageQuanlity\ImageResource\Capture\MotionBlur\'];
img_name = 'Melcaweld01.jpg';
figure('Name', img_name)
for i=1:5
I = imread([img_folder,'L',num2str(i),img_name]);
niqeI = niqe(I,model);
subplot(1,5,i)
imshow(I)
title(['NIQE:',num2str(niqeI)])
end
🌻 计算结果：
3. 使用PIQE
img_folder = ['Z:\LT100\MaMiao\AWorkSpace\项目计划\LT-200\lt-200develop' ...
'\01_可行性研究\OpticalDesign\ImageQuanlity\ImageResource\Capture\GaussianNoise\'];
img_name = 'Melcaweld01.jpg';
figure('Name', img_name)
for i=1:5
I = imread([img_folder,'L',num2str(i),img_name]);
piqeI = piqe(I);
subplot(1,5,i)
imshow(I)
title(['PIQE:',num2str(piqeI)])
end
🌻 计算结果：
4. 结论
- 根据上述实验结果可以知，NIQE算法的评估过程简单，在不同数据的条件下有较好的适应性，评判结果与人主观判断结果较为类似，能够很好地表征图像的质量指标。(作为优选之一)
- 具体的NIQE阈值参数设定需要共同讨论确定，NIQE的值越小，图像的畸变情况越小，目前测试得到LT-200的数据为 9.5 左右，对于 Mecaweld 设备采集的图像其NIQE参数值为 0.9左右（Fitting过程中主要采用了）。(场景的变化存在一定的影响)：
- img_REF0.jpg 图像作为一个优质图像的测试结果也不太理想，但其测试得到的NIQE参数值约为 4.4，因此需要考虑到数据的多样性，还是需要搜素足够的数据量来确保该评估模型的鲁棒性。
- PIQE没有Fitting拟合的过程，通过测试结果可以看出，其对于高斯噪声这一种图像失真的场景得分区分度较差，不如人眼判断的得分情况，不考虑此方法。
🚀 将变换后图像进行混合评分，使用评分机制的图像质量评估算法来进行测试。
二、NIQE 实现
1. 基本代码及注释
calc_NIQE.m
function NIQE = calc_NIQE(input_image_path,shave_width) % 如果不对待进行评估的图像进行 ROI 区域选择，可以设置 shave_width 参数为0
%% Loading model
load modelparameters.mat
blocksizerow = 96; % 切块的行大小
blocksizecol = 96; % 切块的列大小
blockrowoverlap = 0; % 行方向重叠区域大小
blockcoloverlap = 0; % 列方向重叠区域大小
%% Calculating scores
NIQE = [];
input_image = convert_shave_image(imread(input_image_path),shave_width); % 得到经过 [width - 2 * shave_width, height - 2 * shave_width] 裁剪 ROI 之后的灰度图像
% Calculating scores
NIQE = computequality(input_image,blocksizerow,blocksizecol,...
blockrowoverlap,blockcoloverlap,mu_prisparam,cov_prisparam); % 输入待测试图像、分块尺寸、重叠尺寸、mu_prisparam、cov_prisparam，这两个参数的获得从 load modelparameters.mat 加载时就得到了
end
computefeature.m
function feat = computefeature(structdis) % 计算得到的 NSS 模型参数的特征向量
% Input - MSCn coefficients
% Output - Compute the 18 dimensional feature vector
feat = [];
[alpha betal betar] = estimateaggdparam(structdis(:)); % asymmetric generalized Gaussian distribution (AGGD) 计算非对称高斯分布参数
feat = [feat;alpha;(betal+betar)/2];
shifts = [ 0 1;1 0 ;1 1;1 -1];
for itr_shift =1:4 % 在4个偏移方向上对偏移之后的图像进行 AGGD 参数估计
shifted_structdis = circshift(structdis,shifts(itr_shift,:)); % 环形偏移
pair = structdis(:).*shifted_structdis(:); % 计算偏移分布与原分布点积之后的非对称高斯分布参数
[alpha betal betar] = estimateaggdparam(pair); % Eq.7 计算 pair 的非对称高斯分布参数
meanparam = (betar-betal)*(gamma(2/alpha)/gamma(1/alpha)); % Eq.8 分布的均值
feat = [feat;alpha;meanparam;betal;betar];
end
computemean.m
function val = computemean(patch)
val = mean2(patch);
computequality.m
function quality = computequality(im,blocksizerow,blocksizecol,...
blockrowoverlap,blockcoloverlap,mu_prisparam,cov_prisparam)
% Input
% im - Image whose quality needs to be computed
% blocksizerow - Height of the blocks in to which image is divided
% blocksizecol - Width of the blocks in to which image is divided
% blockrowoverlap - Amount of vertical overlap between blocks
% blockcoloverlap - Amount of horizontal overlap between blocks
% mu_prisparam - mean of multivariate Gaussian model 多元高斯模型的均值
% cov_prisparam - covariance of multivariate Gaussian model 多元高斯模型的协方差
% For good performance, it is advisable to use make the multivariate Gaussian model
% using same size patches as the distorted image is divided in to
% Output
%quality - Quality of the input distorted image
% Example call
%quality = computequality(im,96,96,0,0,mu_prisparam,cov_prisparam)
% ---------------------------------------------------------------
%Number of features
% 18 features at each scale
featnum = 18; % 特征空间的特征数量为 18
%----------------------------------------------------------------
%Compute features
if(size(im,3)==3) % 判断图像是否是三通道的数据
im = rgb2gray(im); % 如果是彩色图像则将其转化为灰度图像
end
im = double(im); % 强制转化数据为浮点类型从而方便数据的运算操作，防止操作过程的越界，数据丢失
[row col] = size(im); % 得到灰度图像的 宽 高
block_rownum = floor(row/blocksizerow); % 计算当前图像高度方向上能够分割成为的子图块的数量
block_colnum = floor(col/blocksizecol); % 计算当前图像宽度方向上能够分割成为的子图块的数量
im = im(1:block_rownum*blocksizerow,1:block_colnum*blocksizecol); % 不能分割多余的区域直接舍去，得到新的用于计算的图像区域大小(图像尺寸的自适应配置)
[row col] = size(im); % 得到经过截取选择之后新的图像的大小
block_rownum = floor(row/blocksizerow); % 新的块数量计算
block_colnum = floor(col/blocksizecol);
im = im(1:block_rownum*blocksizerow, ...
1:block_colnum*blocksizecol); % 新的截取之后的图像区域
window = fspecial('gaussian',7,7/6); % 得到一个均值 7 ，方差 7/6 的高斯参数的滤波核
window = window/sum(sum(window)); % 对得到的高斯核进行权重的归一化处理，防止数据越界
scalenum = 2; % 尺寸变化参数设置为 2
warning('off') % 关闭警告
feat = [];
for itr_scale = 1:scalenum
mu = imfilter(im,window,'replicate'); % Eq.2 得到经过高斯滤波之后的图像，这里就得到了 NIQE 算法中的原始图像的 mu均值
mu_sq = mu.*mu;
sigma = sqrt(abs(imfilter(im.*im,window,'replicate') - mu_sq)); % Eq.3 这里得到本地测试图像的 NIQE 算法 sigma标准差，这里对文章中的标准差计算方法进行了简化， 具体的推导可见Readme
structdis = (im-mu)./(sigma+1); % Eq.1 classical spatial NSSmodel 空间 NSS 模型的参数计算，此参数已经能够用来对图像的失真情况进行表征了
feat_scale = blkproc(structdis,[blocksizerow/itr_scale blocksizecol/itr_scale], ...
[blockrowoverlap/itr_scale blockcoloverlap/itr_scale], ...
@computefeature); % 通过对大小为 [m n] 的每个非重叠图像块应用 computefeature 函数并将结果串联到输出矩阵 feat_scale 中，来处理NSS模型参数 structdis。
feat_scale = reshape(feat_scale,[featnum ....
size(feat_scale,1)*size(feat_scale,2)/featnum]);
feat_scale = feat_scale';
if(itr_scale == 1)
sharpness = blkproc(sigma,[blocksizerow blocksizecol], ...
[blockrowoverlap blockcoloverlap],@computemean); % Eq4 计算了所有子块19*19区域的锐度值得平均值，并级联排列在sharpness变量当中
sharpness = sharpness(:);
end
feat = [feat feat_scale];
im =imresize(im,0.5); % 缩放图像到原始图像的一半, 放大时imresize函数使用了双三次插值
end % 每次循环产生有 18 个元素的特征向量，每张子图都需要多计算一次缩放后的特征参数以获得图像的缩放表现，因此共产生 36 个特征值
% Fit a MVG model to distorted patch features
distparam = feat;
mu_distparam = nanmean(distparam); % 用于计算所有去除掉 NaN 之后的参数均值
cov_distparam = nancov(distparam); % 用于计算所有去除掉 NaN 之后的参数协方差
% Compute quality
invcov_param = pinv((cov_prisparam+cov_distparam)/2); % Eq.10 计算矩阵的伪逆
quality = sqrt((mu_prisparam-mu_distparam)* ...
invcov_param*(mu_prisparam-mu_distparam)'); % Eq.10 计算NIQE质量结果
convert_shave_image.m
function shaved = convert_shave_image(input_image,shave_width)
% Converting to y channel only
image_ychannel = rgb2ycbcr(input_image); % 将图像转化为 YCbCr 的数据格式，从而方便提取该格式中的 Y 通道分量
image_ychannel = image_ychannel(:,:,1); % 将 Y 通道分量赋值给 image_ychannel
% Shaving image
shaved = image_ychannel(1+shave_width:end-shave_width,...
1+shave_width:end-shave_width); % 截取一定区域内的灰度图像内容到 shaved 图像当中
end
estimateaggdparam.m
function [alpha betal betar] = estimateaggdparam(vec)
gam = 0.2:0.001:10; % 设置一个从 0.2 变化到 10 的 gamma 参数，步进为 0.001
r_gam = ((gamma(2./gam)).^2)./(gamma(1./gam).*gamma(3./gam)); % gamma 函数是一个扩充到复数域的阶乘计算方法
leftstd = sqrt(mean((vec(vec<0)).^2));
rightstd = sqrt(mean((vec(vec>0)).^2));
gammahat = leftstd/rightstd;
rhat = (mean(abs(vec)))^2/mean((vec).^2);
rhatnorm = (rhat*(gammahat^3 +1)*(gammahat+1))/((gammahat^2 +1)^2);
[min_difference, array_position] = min((r_gam - rhatnorm).^2);
alpha = gam(array_position);
betal = leftstd *sqrt(gamma(1/alpha)/gamma(3/alpha));
betar = rightstd*sqrt(gamma(1/alpha)/gamma(3/alpha));
estimatemodelparam.m
function [mu_prisparam cov_prisparam] = estimatemodelparam(folderpath,...
blocksizerow,blocksizecol,blockrowoverlap,blockcoloverlap,sh_th)
% Input
% folderpath - Folder containing the pristine images
% blocksizerow - Height of the blocks in to which image is divided
% blocksizecol - Width of the blocks in to which image is divided
% blockrowoverlap - Amount of vertical overlap between blocks
% blockcoloverlap - Amount of horizontal overlap between blocks
% sh_th - The sharpness threshold level
%Output
%mu_prisparam - mean of multivariate Gaussian model
%cov_prisparam - covariance of multivariate Gaussian model
% Example call
%[mu_prisparam cov_prisparam] = estimatemodelparam('pristine',96,96,0,0,0.75);
%----------------------------------------------------------------
% Find the names of images in the folder
current = pwd;
cd(sprintf('%s',folderpath))
names = ls;
names = names(3:end,:);
cd(current)
% ---------------------------------------------------------------
%Number of features
% 18 features at each scale
featnum = 18;
% ---------------------------------------------------------------
% Make the directory for storing the features
mkdir(sprintf('local_risquee_prisfeatures'))
% ---------------------------------------------------------------
% Compute pristine image features
for itr = 1:size(names,1)
itr
im = imread(sprintf('%s\\%s',folderpath,names(itr,:)));
if(size(im,3)==3)
im = rgb2gray(im);
end
im = double(im);
[row col] = size(im);
block_rownum = floor(row/blocksizerow);
block_colnum = floor(col/blocksizecol);
im = im(1:block_rownum*blocksizerow, ...
1:block_colnum*blocksizecol);
window = fspecial('gaussian',7,7/6);
window = window/sum(sum(window));
scalenum = 2;
warning('off')
feat = [];
for itr_scale = 1:scalenum
mu = imfilter(im,window,'replicate');
mu_sq = mu.*mu;
sigma = sqrt(abs(imfilter(im.*im,window,'replicate') - mu_sq));
structdis = (im-mu)./(sigma+1);
feat_scale = blkproc(structdis,[blocksizerow/itr_scale blocksizecol/itr_scale], ...
[blockrowoverlap/itr_scale blockcoloverlap/itr_scale], ...
@computefeature);
feat_scale = reshape(feat_scale,[featnum ....
size(feat_scale,1)*size(feat_scale,2)/featnum]);
feat_scale = feat_scale';
if(itr_scale == 1)
sharpness = blkproc(sigma,[blocksizerow blocksizecol], ...
[blockrowoverlap blockcoloverlap],@computemean);
sharpness = sharpness(:);
end
feat = [feat feat_scale];
im =imresize(im,0.5);
end
save(sprintf('local_risquee_prisfeatures\\prisfeatures_local%d.mat',...
itr),'feat','sharpness');
end
%----------------------------------------------
% Load pristine image features
prisparam = [];
current = pwd;
cd(sprintf('%s','local_risquee_prisfeatures'))
names = ls;
names = names(3:end,:);
cd(current)
for itr = 1:size(names,1)
% Load the features and select the only features
load(sprintf('local_risquee_prisfeatures\\%s',strtrim(names(itr,:))));
IX = find(sharpness(:) >sh_th*max(sharpness(:)));
feat = feat(IX,:);
prisparam = [prisparam; feat];
end
%----------------------------------------------
% Compute model parameters
mu_prisparam = nanmean(prisparam);
cov_prisparam = nancov(prisparam);
%----------------------------------------------
% Save features in the mat file
save('modelparameters_new.mat','mu_prisparam','cov_prisparam');
%----------------------------------------------
2. 参考文献
- 《Making a (Completely Blind) Image Quality Analyzer.pdf》
- 《Multiscale_skewed_heavy_tailed_model_for_texture_analysis.pdf》
三、基础数学原理
1、广义高斯分布(GGD)
1.1 描述
零均值的广义高斯分布如下：
\[f(x;\gamma,\sigma^2)=\frac{\gamma}{2\beta\Gamma(1/\gamma)}exp[-(\frac{|x|}{\beta})^\gamma]
\]
其中\(\beta=\sigma\sqrt\frac{\Gamma(1/\gamma)}{\Gamma(3/\gamma)}\)
而 \(\Gamma(x)\) 称为 Gamma 函数，其递归应用满足 \(\Gamma(n-1)=\frac{\Gamma(n)}{n-1}\) ，形状参数 \(\gamma\) 控制了分布的形状，而 \(\sigma^2\) 控制了方差，例如取 \(\gamma = 2\) 就会得到零均值的高斯分布：
\[\beta=\sigma\sqrt\frac{\Gamma(1/2)}{\Gamma(3/2)}=\sqrt2 \sigma
\]
\[\begin{split}
f(x)
&=\frac{\gamma}{2\beta\Gamma(1/\gamma)}exp[-(\frac{|x|}{\beta})^\gamma] \\
&=\frac{1}{\sqrt2\sigma\sqrt\pi}exp[-(\frac{|x|}{\sqrt2\sigma})^2] \\
&=\frac{1}{\sqrt{2\pi}\sigma}exp[-\frac{x^2}{2\sigma^2}]
\end{split}
\]
首先记
\(a=\frac{\gamma}{2\beta\Gamma(1/\gamma)}, b=\frac{1}{\beta}\)
则
\[E[|X|]=\int^{+\infin}_{-\infin}|x|f_X(x)dx=\int^{+\infin}_{-\infin}|x|ae^{-[b|x|]^\gamma}dx
\]
取 \(g(x) = |x|ae^{-[b|x|]^\gamma} = g(-x)\)，即 \(g(x)\) 为偶函数，正半轴和负半轴上的积分结果相等，上式简化如下：
\[\begin{split}
E[|X|]
&=\int^{+\infin}_{0}2axe^{-(bx)^\gamma}dx, (bx)^\gamma = y \\
&=\frac{2a}{b^2\gamma}\int^{+\infin}_{0}\gamma^{(2/\gamma)-1}e^{-y}dy=\frac{2a}{b^2\gamma}\Gamma(2/\gamma)
\end{split}
\]
因此
\[E[|X|]=\frac{2a}{b^2\gamma}\Gamma(2/\gamma)=\frac{\beta}{\Gamma(1/\gamma)}=\sigma\sqrt\frac{\Gamma(1/\gamma)}{\Gamma(3/\gamma)}\frac{\Gamma(2/\gamma)}{\Gamma(1/\gamma)}
\]
则
\[E^2[|X|]=\sigma^2\frac{\Gamma^2(2/\gamma)}{\Gamma(3/\gamma)\Gamma(1/\gamma)}
\]
得到一个比函数：
\[\rho(\gamma)=\frac{\sigma^2}{E^2[|X|]}=\frac{\Gamma(3/\gamma)\Gamma(1/\gamma)}{\Gamma^2(2/\gamma)}
\]
1.2 参数估计方法
对于零均值广义高斯分布，计算估计值：
\[\hat{E}[|X|]=\frac{1}{M}\sum^M_{i=1}|x_i| \\
\hat{\sigma}^2_X=\frac{1}{M}\sum^M_{i=1}x_i^2
\]
根据上面的结果计算可得：
\[\frac{\hat{\sigma}^2_X}{\hat{E}[|X|]}=\hat{\rho}(\gamma)=\frac{\Gamma(3/\hat{\gamma})\Gamma(1/\hat{\gamma)}}{\Gamma^2(2/\hat{\gamma)}}
\]
直到了 \(\rho(\gamma)\) 的估计值之后，就很容易通过枚举方式来估计 \(\gamma\)
1.3 代码
function [gamparam sigma] = estimateggdparam(vec)
gam = 0.2:0.001:10;
r_gam = (gamma(1./gam).*gamma(3./gam))./((gamma(2./gam)).^2);
sigma_sq = mean((vec.^2));
sigma = sqrt(sigma_sq);
E = mean(abs(vec));
rho = sigma_sq/E^2;
[min_difference, array_position] = min(abs(rho - r_gam));
gamparam = gam(array_position);
2、非对称广义高斯分布(ADDG)
2.1 描述
零均值的非对称广义高斯分布如下：
\[f(x;\alpha,\sigma^2_l,\sigma^2_r)=
\begin{cases}
\frac{\alpha}{(\beta_l+\beta_r)\Gamma(1/\alpha)}exp(-(\frac{-x}{\beta_l})^\alpha)\quad x<0 \\
\frac{\alpha}{(\beta_l+\beta_r)\Gamma(1/\alpha)}exp(-(\frac{x}{\beta_r})^\alpha)\quad x>0
\end{cases}
\]
其中
\[\beta_l = \sigma_l\sqrt{\frac{\Gamma(1/\alpha)}{\Gamma(3/\alpha)}},\quad \beta_r = \sigma_r\sqrt{\frac{\Gamma(1/\alpha)}{\Gamma(3/\alpha)}}
\]
形状参数 \(\alpha\) 控制分布的形状，而 \(\sigma^2_l\) 和 \(\sigma^2_r\) 是缩放参数，他们控制模式两边的扩散程度，当 \(\sigma^2_l=\sigma^2_r\) 时，AGGD 退化成 GGD，这里可以参考 《Multiscale Skewed Heavy Tailed Model For Texture Analysis》 .
记
\[y=\mu x^\gamma
\]
则
\[x=(y/\mu)^{1/\gamma},dx=\frac{1}{\gamma}(y/\mu)^{(1/y)-1}\frac{1}{\mu}dy
\]
因此
\[\begin{split}
\int^{+\infin}_0x^{v-1}exp(-\mu x^\gamma)dx
&= \int^{+\infin}_0(y/\mu)^{(v-1)/\gamma}e^{-y}\frac{1}{\gamma}(y/\mu)^{(1/y)-1}\frac{1}{\mu}dy\\
&= \frac{1}{\gamma \mu}\int^{+\infin}_0(y/\mu)^{v/\gamma -1}e^{-y}dy \\
&= \frac{1}{\gamma \mu}\mu^{1-v/\gamma}\int^{+\infin}_0y^{v/\gamma-1}e^{-y}dy \\
&= \frac{1}{\gamma}\mu^{-v/\gamma}\cdot \Gamma(v/\gamma)
\end{split}
\]
记
\[a=\frac{\alpha}{(\beta_l+\beta_r)\Gamma(1/\alpha)}, b=\frac{1}{\beta_l}, c=\frac{1}{\beta_r}
\]
就有
\[\begin{split}
\mu_k = E[X^k]=\int^{+\infin}_{-\infin}x^kf_X(x)dx
&= \int^{0}_{-\infin}x^kf_X(x)dx+ \int^{+\infin}_0x^kf_X(x)dx\\
&= \int^{0}_{-\infin}x^kaexp(-(-bx)^\alpha)dx+ \int^{+\infin}_0x^kaexp(-(cx)^\alpha)dx\\
&= (-1)^{k+1}a\int^{0}_{+\infin}x^kexp(-(bx)^\alpha)dx+ \int^{+\infin}_0x^kaexp(-(cx)^\alpha)dx \quad \text{use substitution method:} x=-t\\
&= (-1)^{k}a\int^{+\infin}_{0}x^kexp(-(bx)^\alpha)dx+ \int^{+\infin}_0x^kaexp(-(cx)^\alpha)dx\\
&= (-1)^k\frac{a}{\alpha}b^{-(k+1)}\cdot \Gamma(\frac{k+1}{\alpha}) + \frac{a}{\alpha}c^{-(k+1)}\cdot \Gamma(\frac{k+1}{\alpha}) \\
&= \frac{a}{\alpha}\cdot \Gamma(\frac{k+1}{\alpha})[(-1)^kb^{-(k+1)} + c^{-(k+1)}]=\frac{a}{\alpha}\cdot \Gamma(\frac{k+1}{\alpha})[(-1)^k\beta_l^{k+1} + \beta_r^{k+1}]\\
&= \frac{(-1)^k\cdot \beta_l^{k+1} + \beta_r^{k+1}}{\beta_l+\beta_r}\frac{\Gamma((k+1)/\alpha)}{\Gamma(1/\alpha)}
\end{split}
\]
根据上述推导结果可知：
\[\begin{split}
m_k = E[|X|^k]=\int^{+\infin}_{-\infin}|x|^kf_X(x)dx
&= \frac{\beta_l^{k+1} + \beta_r^{k+1}}{\beta_l+\beta_r}\frac{\Gamma((k+1)/\alpha)}{\Gamma(1/\alpha)}
\end{split}
\]
然后计算比值：
\[\begin{split}
r=\frac{m_1^2}{\mu_2}
&=\frac{(\beta^2_l+\beta^2_r)^2}{(\beta_l+\beta_r)}\frac{\Gamma^2(2/\alpha)}{\Gamma^2(1/\alpha)}\frac{\beta_l+\beta_r}{\beta^3_l+\beta^3_r}\frac{\Gamma(1/\alpha)}{\Gamma(3/\alpha)}\\
&= \frac{\Gamma^2(2/\alpha)}{\Gamma(1/\alpha)\Gamma(3/\alpha)}\frac{(\beta^2_l+\beta^2_r)^2}{(\beta_l+\beta_r)(\beta^3_l+\beta^3_r)}\\
&= \frac{(\frac{\beta^2_l+\beta^2_r}{\beta^2_r})^2}{(\frac{\beta_l+\beta_r}{\beta_r})(\frac{\beta^3_l+\beta^3_r}{\beta^3_r})}\frac{\Gamma^2(2/\alpha)}{\Gamma(1/\alpha)\Gamma(3/\alpha)}=\frac{(\gamma^2+1)^2}{(\gamma+1)(\gamma^3+1)}\cdot \rho(\alpha)
\end{split}
\]
其中 \(\gamma=\frac{\beta_l}{\beta_r}, \rho(\alpha)=\frac{\Gamma^2(2/\alpha)}{\Gamma(1/\alpha)\Gamma(3/\alpha)}\)
2.2 参数估计方法
首先估计 \(\sigma_r\) 和 \(\sigma_l\) ：
\[\hat{\sigma_l^2}=\frac{1}{N_l-1}\sum^{N_l}_{k=1,x_k<0}x_k^2, \quad \hat{\sigma_r^2}=\frac{1}{N_r-1}\sum^{N_r}_{k=1,x_r<0}x_k^2
\]
所以
\[\gamma=\frac{\beta_l}{\beta_r}=\frac{\sigma_l}{\sigma_r}, \quad \hat{\gamma}=\frac{\hat{\sigma_l}}{\hat{\sigma_r}}=\sqrt{\frac{\frac{1}{N_l-1}\sum^{N_l}_{k=1,x_k<0}x_k^2}{\frac{1}{N_r-1}\sum^{N_r}_{k=1,x_r<0}x_k^2}}
\]
而 \(r\) 的一个无偏估计是
\[\hat{r}=\frac{\hat{m}_1^2}{\hat{\mu}_2}=\frac{(\sum|x_k|)^2}{\sum x_k^2}
\]
所以就可以得到
\[r=\frac{(\gamma^2+1)^2}{(\gamma+1)(\gamma^3+1)}\cdot\rho(\alpha)
\]
求得
\[\hat{\rho}(\alpha)=\frac{(\hat{\gamma}+1)(\hat{\gamma}^3+1)}{(\hat{\gamma}^2+1)^2}\cdot \hat{r}
\]
2.3 代码
来自 BRISQUE 的 Matlab 代码
function [alpha leftstd rightstd] = estimateaggdparam(vec)
gam = 0.2:0.001:10;
r_gam = (gamma(2./gam)).^2./(gamma(1./gam).*gamma(3./gam));
leftstd = sqrt(mean((vec(vec<0)).^2));
rightstd = sqrt(mean((vec(vec>0)).^2));
gammahat = leftstd / rightstd;
rhat = (mean(abs(vec)))^2/mean((vec).^2);
rhatnorm = (rhat*gammahat^3 + 1) * (gammahat + 1) / ((gammahat^2 + 1)^2);
[min_difference, array_position] = min((r_gam - rhatnorm).^2);
alpha = gam(array_position);
Reference
- 参考博客：https://blog.csdn.net/qq_30124657/article/details/114268329
- python代码：https://github.com/guptapraful/niqe
- 《Learning a No-Reference Quality Metric for Single-Image Super-Rolution", CVIU 2017.》
- matlab的代码，路径是：https://github.com/chaoma99/sr-metric
- NIQE推导
- Multiscale_skewed_heavy_tailed_model_for_texture_analysis.pdf
转载请注明出处！感谢GISPALAB实验室的老师和同学们的帮助和支持~
浙公网安备 33010602011771号