圖像質量評估-BRISQUE

1. 介紹

1.1 什么是圖像質量評估Image Quality Assessment (IQA)？

圖像質量評價(IQA)算法以任意圖像作為輸入，輸出質量分數作為輸出。有三種類型的IQA：

• 全參考圖像質量評價
  適用情形：您有一個“干凈”參考(非扭曲)圖像，以衡量您的扭曲圖像的質量。此度量可用於評估圖像壓縮算法的質量，在該算法中，我們可以訪問原始圖像及其壓縮版本。
• 半參考圖像質量評價
  適用情形：如果沒有參考圖像，而是具有一些選擇性信息的圖像(例如，水印圖像)來比較和測量失真圖像的質量。
• 無參考圖像質量評價
  適用情形：算法得到的唯一輸入是要測量其質量的圖像

1.2 無參考圖像質量評價

在這篇文章中，我們將討論一個無參考圖像質量評價的IQA度量算法，稱為盲/無參考圖像空間質量評估器Blind/Referenceless Image Spatial Quality Evaluator (BRISQUE)。在深入研究這個理論之前，讓我們先了解兩個基本術語：

• 失真圖像(扭曲圖像)
  顧名思義，失真圖像是指被模糊、噪聲、水印、顏色變換、幾何變換等因素扭曲的原始圖像。下圖是TID 2008數據庫中使用的圖像失真情況。
  
• 自然圖像
  自然圖像是指由相機直接拍攝沒有后期處理的圖像。以下自然圖像和失真圖像的示例。左圖是自然圖像，右圖是失真圖像。
  
  正如你所能想象的，圖像清晰度與圖像是失真的還是自然的沒有關系。例如，當視頻被巧妙地用運動模糊渲染時，由於故意模糊，算法可能會對其質量產生混淆。因此，我們必須在正確的背景下使用這種正確的圖像質量評價方法來度量圖像。

1.3 圖像質量評估(IQA)數據集

圖像質量好於壞是一個主觀問題。為了獲得一個優秀的圖像質量評估算法，我們需要給出許多圖像的算法示例和它們的質量分數。誰為這些訓練圖像給定質量得分？當然是人類。但我們不能僅僅依靠一個人的意見。因此，我們需要綜合評估個人的意見，並為圖像分配0（最佳）和100（最差）之間的平均分數。該分數在學術文獻中稱為平均質量分數。幸運的是我們不需要自己收集這些數據，有一個數據集稱為TID 2008已經被用於研究目的。有關TID2008的信息見：
https://pdfs.semanticscholar.org/6dd3/7b57f3391b438fa588f98a1c2067365ae5ca.pdf

如下圖所示，TID2008圖像質量得分范圍為0到100，得分越小，主觀質量越好。

2. 盲/無參考圖像空間質量評估器（BRISQUE）

下圖給出了計算BRISQUE所涉及的步驟

主要為三步：

1. 提取自然場景統計（NSS）
2. 計算特征向量
3. 預測圖像質量得分

2.1. 提取自然場景統計（NSS）

自然圖像的像素強度的分布不同於失真圖像的像素強度的分布。當歸一化像素強度並計算這些歸一化強度上的分布時，這種分布差異更加明顯。特別是經過歸一化后的自然圖像像素強度服從高斯分布(Bell曲線)，而非自然圖像或畸變圖像的像素強度不服從高斯分布。因此，通過像素高斯分布曲線是衡量圖像失真的一種方法。我們已在下圖中說明了這一點。

上圖左圖顯示了一幅有添加人工效應的自然圖像，符合高斯分布。右邊為一個人造圖像，但不適合同樣的分布。

2.1.1 MSCN（Mean Subtracted Contrast Normalized）系數

有幾種不同的方法來規范化圖像。一種這樣的歸一化稱為MSCN。下圖顯示了如何計算MSCN系數。

上圖MSCN計算過程能夠可視化為：

下面要講述詳細的數學推導，但不要讓下面的數學嚇倒你。數學之后的代碼要簡單得多，容易理解！

代碼如下：
C++：

Mat im = imread("image_scenery.jpg"); // read image
cvtColor(im, im, COLOR_BGR2GRAY); // convert to grayscale

im.convertTo(im, 1.0/255); // normalize and copy the read image to orig_bw
Mat mu(im.size(), CV_64FC1, 1);
GaussianBlur(im, mu, Size(7, 7), 1.166); // apply gaussian blur
Mat mu_sq = mu.mul(mu); 

// compute sigma
Mat sigma = im.size(), CV_64FC1, 1);
sigma = im.mul(im); 
GaussianBlur(sigma, sigma, Size(7, 7), 1.166); // apply gaussian blur
subtract(sigma, mu_sq, sigma); // sigma = sigma - mu_sq
cv::pow(sigma, double(0.5), sigma); // sigma = sqrt(sigma)
add(sigma, Scalar(1.0/255), sigma); // to avoid DivideByZero Exception
Mat structdis(im.size(), CV_64FC1, 1);
subtract(im, mu, structdis); // structdis = im - mu
divide(structdis, sigma, structdis); 

Python：

im = cv2.imread("image_scenery.jpg", 0) # read as gray scale
blurred = cv2.GaussianBlur(im, (7, 7), 1.166) # apply gaussian blur to the image
blurred_sq = blurred * blurred 
sigma = cv2.GaussianBlur(im * im, (7, 7), 1.166) 
sigma = (sigma - blurred_sq) ** 0.5
sigma = sigma + 1.0/255 # to make sure the denominator doesn't give DivideByZero Exception
structdis = (im - blurred)/sigma # final MSCN(i, j) image

2.1.2 相鄰像素間的乘積關系

MSCN為像素強度提供了不錯的歸一化。然而自然退休與失真圖像的差異不僅限於像素強度分布，還包括相鄰像素之間的關系。為了捕獲相鄰像素的關系，在四個方向來求出相鄰元素的兩兩乘積，即：水平（H），垂直（V），左對角線（D1），右對角線（D2）。

可以用Python和C++計算兩兩乘積，如下所示：
C++：

// declare shifting indices array
int shifts[4][2] = {{0, 1}, {1, 0}, {1, 1}, {-1, 1}};

// calculate pair-wise products for every combination of shifting indices
for(int itr_shift = 1; itr_shift <= 4; itr_shift++) 
{
    int* reqshift = shifts[itr_shift - 1]; // the required shift index
    // declare shifted image
    Mat shifted_structdis(imdist_scaled.size(), CV_64F, 1);
    // BwImage is a helper class to create a copy of the image and create helper functions to access it's pixel values
    BwImage OrigArr(structdis);
    BwImage ShiftArr(shifted_structdis); 
    // calculate pair-wise component for the given orientation
    for(int i = 0; i < structdis.rows; i++)
    {
        for(int j = 0; j < structdis.cols; j++) { if(i + reqshift[0] >= 0 && i + reqshift[0] < structdis.rows && j + reqshift[1] >= 0 && j + reqshift[1] < structdis.cols)
            {
                ShiftArr[i][j] = OrigArr[i + reqshift[0]][j + reqshift[1]];
            }f
            else
            {
                ShiftArr[i][j] = 0;
            }
        }   
    }
    Mat shifted_new_structdis;
    shifted_new_structdis = ShiftArr.equate(shifted_new_structdis);
    // find the pairwise product
    multiply(structdis, shifted_new_structdis, shifted_new_structdis);
}

Python：

# indices to calculate pair-wise products (H, V, D1, D2)
shifts = [[0,1], [1,0], [1,1], [-1,1]]
# calculate pairwise components in each orientation
for itr_shift in range(1, len(shifts) + 1):
    OrigArr = structdis
    reqshift = shifts[itr_shift-1] # shifting index

    for i in range(structdis.shape[0]):
        for j in range(structdis.shape[1]):
            if(i + reqshift[0] >= 0 and i + reqshift[0] < structdis.shape[0] \ and j + reqshift[1] >= 0 and j  + reqshift[1] < structdis.shape[1]):
               ShiftArr[i, j] = OrigArr[i + reqshift[0], j + reqshift[1]]
            else:
               ShiftArr[i, j] = 0

可以使用cv2.warpAffine方法將兩個for循環簡化為幾行代碼。這將大大加快計算速度。

# create affine matrix (to shift the image)
M = np.float32([[1, 0, reqshift[1]], [0, 1, reqshift[0]]])
ShiftArr = cv2.warpAffine(OrigArr, M, (structdis.shape[1], structdis.shape[0])

2.2 計算特征向量

到目前為止，我們已經從原始圖像派生了5個圖像參數，即1個MSCN圖像和4個兩兩乘積圖像，以捕獲相鄰關系(水平、垂直、左對角、右對角)。接下來，我們將使用以上5個參數來計算尺寸為36×1的特征向量。請注意，原始輸入圖像可以是任何尺寸（寬度/高度），但特征向量的大小始終為36×1。
通過將MSCN圖像擬合到廣義高斯分布(GGD)來計算36×1特征向量的前兩個元素。GGD有兩個參數-一個用於形狀，另一個用於方差。
接下來，用非對稱廣義高斯分布（AGGD）擬合相鄰元素乘積參數中的任何一個。AGGD是廣義高斯擬合（GGD）的不對稱形式。它有四個參數，即形狀，均值，左方差和右方差。由於有4個兩兩相乘圖像參數，我們最終得到16個值。
因此，我們最終得到了特征向量的18個元素。圖像縮小到原來大小的一半，並重復相同的過程以獲得18個新的數字，使總數達到36個。下表概述了這一點。

2.3 圖像質量評分的預測

在典型的機器學習應用程序中，首先將圖像轉換為特征向量。然后將訓練數據集中所有圖像的特征向量和輸出(在這種情況下是質量分數)輸入到支持向量機(SVM)等學習算法中。可以下載tid2008數據並且訓練一個支持向量機來解決這個問題，但是在這篇文章中，我們將簡單地使用作者提供的經過訓練的模型。
下載地址：http://www.ponomarenko.info/tid2008.htm
首先加載訓練后的模型，然后利用模型產生的支持向量預測概率，利用LIBSVM預測最終的質量分數。需要注意的是，特征向量首先被縮放為-1到1，然后用於預測。我們共享Python和C++實現這一點的方法：

3. 結果和代碼

3.1 結果

我們對四種類型的失真進行了度量。以下是理論上每個失真的最終質量分數。但是這個僅僅是理論上的，實際上效果在不同電腦和不同圖像分辨率上並不是這樣的。對於圖像壓縮結果可能大相徑庭。不過總的來說BRISQUE效果不錯，但是精度不夠，主要是svm還需要自己訓練。

3.2 代碼

TODO

4. 最后