MSCN（Mean Subtracted Contrast Normalized）系數的直方圖

      MSCN系數是無參考的空間域圖像質量評估算法BRISQUE（No-Reference Image Quality Assessment in the Spatial Domain）中提出的，MSCN系數具有由於失真的存在而改變的特征統計特性，並且量化這些變化將使得可以預測影響圖像的失真類型以及其感知質量。這篇論文的大致原理是從圖像中提取MSCN系數，然后將MSCN系數擬合成非對稱性廣義高斯分布，提取擬合的高斯分布的特征，輸入到支持向量機SVM中做回歸，最終得到圖像質量的評分。

       這篇論文提供了源代碼，其中有MSCN的計算代碼，但是論文中的關於證明不同失真圖像的MSCN系數具有不同分布的直方圖的代碼並沒有，經過大概大半天的時間，我把這個直方圖用matlab做了出來。

      首先，直觀的看一下MSCN系數圖如下（原圖和MSCN系數圖）：

       

                                                                                                                                                

MSCN系數圖的計算公式源論文中有，代碼從源代碼中修改而來，如下：

clear;clc;
imdist=imread('2.png');

%imshow(imdist); %顯示原圖
imdist = rgb2gray(imdist);
imdist = double(imdist);

%imdist=phasecong3(imdist);
%imshow(phaseCong)

window = fspecial('gaussian',7,7/6);
window = window/sum(sum(window));

mu            = filter2(window, imdist, 'same');
mu_sq         = mu.*mu;
sigma         = sqrt(abs(filter2(window, imdist.*imdist, 'same') - mu_sq));
structdis     = (imdist-mu)./(sigma+1);


%imshow(structdis,[]);  %顯示MSCN圖

　　其實MSCN系數計算出來是一個浮點數矩陣，就是上面代碼中的structdis，范圍是-2到2。論文中的SMCN系數直方圖橫坐標表示圖像的MSCN系數值，縱坐標表示值出現的頻率，在網上找了好久，發現並沒有源代碼，只好自己寫了。剛開始的時候以為每個值會出現多次，統計了一下發現絕大多數MSCN值僅出現了一次，后來才想明白，橫坐標表示的不是每一個MSCN值的頻率，而是一個范圍值的頻率。浮點數不好統計，於是剛開始的時候我把MSCN系數的值轉化到了0-255，然后再來統計它的頻率，統計個數用的是matlab自帶的函數tabulate，做出的圖如下：

具體代碼如下：

clear;clc;
imdist=imread('2.png');

%imshow(imdist);
imdist = rgb2gray(imdist);
imdist = double(imdist);

%imdist=phasecong3(imdist);
%imshow(phaseCong)

window = fspecial('gaussian',7,7/6);
window = window/sum(sum(window));

mu            = filter2(window, imdist, 'same');
mu_sq         = mu.*mu;
sigma         = sqrt(abs(filter2(window, imdist.*imdist, 'same') - mu_sq));
structdis     = (imdist-mu)./(sigma+1);


%imshow(structdis,[]);  %顯示MSCN圖
% imshow(mu,[]);
% imshow(sigma,[]);
%hist(structdis);

  InImg=structdis;
  ymax=255;ymin=0;
  xmax = max(max(InImg)); %求得InImg中的最大值
  xmin = min(min(InImg)); %求得InImg中的最小值
  OutImg = round((ymax-ymin)*(InImg-xmin)/(xmax-xmin) + ymin); %歸一化並取整

%hist(structdis);

n=tabulate(OutImg(:)); %計算矩陣中各元素出現的頻數
number=n(:,2);

a=number;

x=1:1:245;
plot(x,a,'-*b'); %繪圖

　　其他問題先不說，首先就是點太密集了，跟論文中的完全不一樣，而且MSCN值轉化到0-255之后可能會丟失一些精度。所以這個不合格。然后原因出在matlab自帶的統計矩陣值個數的函數上，這個函數不靠譜，只能自己寫了，倒不是很麻煩，循環統計范圍值就好了，但是以前寫c++習慣了，用matlab寫的時候總是想着c++的方式來寫，結果很簡單的問題頭都搞暈了，還好有好朋友幫忙，循環很快寫出來了。然后封裝成函數如下（以下的圖來自TID2013的bike.bmp）：

function [n]=compute_mscn(imdist)    %計算圖像的MSCN系數並統計其像素值出現的頻數

%imdist=imread('3.png');
imdist = rgb2gray(imdist);
imdist = double(imdist);

%imdist=phasecong3(imdist);
%imshow(phaseCong)

window = fspecial('gaussian',7,7/6);
window = window/sum(sum(window));

mu            = filter2(window, imdist, 'same');
mu_sq         = mu.*mu;
sigma         = sqrt(abs(filter2(window, imdist.*imdist, 'same') - mu_sq));
structdis     = (imdist-mu)./(sigma+1);
% imshow(structdis,[]);  %顯示MSCN圖
% imshow(mu,[]);
% imshow(sigma,[]);
%hist(structdis);

  InImg=structdis;
  ymax=255;ymin=0;
  xmax = max(max(InImg)); %求得InImg中的最大值
  xmin = min(min(InImg)); %求得InImg中的最小值
  OutImg = round((ymax-ymin)*(InImg-xmin)/(xmax-xmin) + ymin); %歸一化並取整

%hist(structdis);

%n=tabulate(OutImg(:)); %計算矩陣中各元素出現的頻數
%number=n(:,2);

%a=number;
for j=1:52         %統計MSCN系數出現的頻數
    i=(j-1)*5;
    n(j)=length(find(OutImg>=i&OutImg<i+2));
    % t=t+1;
end
end

　　調用的代碼也很簡單，如下所示：

imdist=imread('bikes.bmp');
n1=compute_mscn(imdist);

imdist=imread('img15.bmp');
n2=compute_mscn(imdist);

imdist=imread('img31.bmp');
n3=compute_mscn(imdist);

imdist=imread('img129.bmp');
n4=compute_mscn(imdist);

x=0:5:255;
%x=-2:0.05:2;
plot(x,n1,'-*b',x,n2,'-dr',x,n3,'-sg',x,n4,'-ok');
set(gca,'XTick',[0:50:250])
%set(gca,'XTick',[-2:0.5:2]);%x軸范圍，間隔0.5
legend('org','jp2k','ff','gblur');   %右上角標注
xlabel('MSCN系數')  %x軸坐標描述
ylabel('頻數（歸一化）') %y軸坐標描述
grid on;  %加網格線
%set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1); %加網格虛線

　　然后結果圖如下：

 

理論上和論文中的意思是一樣了，但是還是有些差距，首先是MSCN系數的值被歸一化到0-255了，然后頻數還沒有歸一化，所以考慮不把MSCN系數的值轉換到0-255，直接進行統計，以免有誤差，代碼仍然進行了封裝，函數如下：

function [n]=compute_mscn2(imdist)    %計算圖像的MSCN系數並統計其像素值出現的頻數

imdist = rgb2gray(imdist);
imdist = double(imdist);

%imdist=phasecong3(imdist);
%imshow(phaseCong)

window = fspecial('gaussian',7,7/6);
window = window/sum(sum(window));

mu            = filter2(window, imdist, 'same');
mu_sq         = mu.*mu;
sigma         = sqrt(abs(filter2(window, imdist.*imdist, 'same') - mu_sq));
structdis     = (imdist-mu)./(sigma+1);
% imshow(structdis,[]);  %顯示MSCN圖
% imshow(mu,[]);
% imshow(sigma,[]);
%hist(structdis);

for j=1:81              %統計MSCN系數出現的頻數
    i=(j-40)*0.05;
    n(j)=length(find(structdis >=i&structdis<i+0.05));
    % t=t+1;
end
n=mapminmax(n,0,1);

end

　　調用的代碼如下：

imdist=imread('bikes.bmp');
n1=compute_mscn2(imdist);

imdist=imread('img15.bmp');
n2=compute_mscn2(imdist);

imdist=imread('img31.bmp');
n3=compute_mscn2(imdist);

imdist=imread('img129.bmp');
n4=compute_mscn2(imdist);

%x=0:5:255;
x=-2:0.05:2;
plot(x,n1,'-*b',x,n2,'-dr',x,n3,'-sg',x,n4,'-ok');
%set(gca,'XTick',[0:50:250])
set(gca,'XTick',[-2:0.5:2]);%x軸范圍，間隔0.5
legend('org','jp2k','ff','gblur');   %右上角標注
xlabel('MSCN系數')  %x軸坐標描述
ylabel('頻數（歸一化）') %y軸坐標描述
grid on;  %加網格線
%set(gca,'GridLineStyle',':','GridColor','k','GridAlpha',1); %加網格虛線

　　最終總算是完成任務了，可以看出確實不同失真圖像的MSCN系數具有不同的直方圖分布，圖如下：