一、簡介

方向梯度直方圖(Histogram of Oriented Gradient,HOG)是用於在計算機視覺和影象處理領域，目標檢測的特徵描述子。該項技術是用來計算影象區域性出現的方向梯度次數或資訊進行計數。此種方法跟邊緣方向直方圖、尺度不變特徵變換以及形狀上下文方法有很多相似。但與它們的不同點是：HOG的計算基於一致空間的密度矩陣來提高準確率。即：在一個網格密集的大小統一的細胞單元上計算，而且為了提高效能，還採用了重疊的區域性對比度歸一化技術。HoG特徵與SVM分類器結合，已經被廣泛應用於影象識別中，尤其在行人檢測。本節內容介紹HOG相關知識。

1 基本介紹

提出HOG是由Navneet Dalal & Bill Triggs在2005年發表在CVPR中論文[1]。作者給了特徵提取的流程圖，如下圖1所示：

HOG的核心思想是所檢測的區域性物體外形能夠被光強梯度或邊緣方向的分佈所描述。通過將整幅影象分割成小的連線區域稱為cells，每個cell生成一個方向梯度直方圖或者cell中pixel的邊緣方向，這些直方圖的組合可表示出所檢測目標的目標）描述子。為改善準確率，區域性直方圖可以通過計算影象中一個較大區域稱為block的光強作為measure被對比標準化，然後用這個measure歸一化這個block中的所有cells.這個歸一化過程完成了更好的照射/陰影不變性。與其他描述子相比，HOG得到的描述子保持了幾何和光學轉化不變性除非物體方向改變。而Block與Cells關係如下圖所示：

而現在，我們給出作者做的行人檢測試驗，如下圖6所示：

其中，圖中(a)表示所有訓練影象集的平均梯度；(b)和(c)分別表示：影象中每一個區間上的最大最大正、負SVM權值；(d)表示一副測試影象；(e)計算完R-HOG後的測試影象；(f)和(g)分別表示被正、負SVM權值加權後的R-HOG影象。

3 演算法描述

HOG特徵提取方法：首先將影象分成小的連通區域，我們把它叫細胞單元。然後採集細胞單元中各畫素點的梯度的或邊緣的方向直方圖。最後把這些直方圖組合起來就可以構成特徵描述器。如下圖所示。

即：將一個影象：第一步，灰度化即：將影象看做一個x,y,z(灰度)的三維影象；第二步：劃分成小cells(2*2)；第三步，計算每個cell中每個pixel的gradient(即orientation)；最後， 統計每個cell的梯度直方圖(不同梯度的個數)，即可形成每個cell的descriptor。整個演算法具體過程由以下幾個部分組成：

4 色彩和伽馬歸一化(Gamma/Color normalization)

在實際在實際應用中，這一步可以省略，因為作者分別在灰度空間、RGB色彩空間和LAB色彩空間上對影象進行色彩和伽馬歸一化。而實驗結果顯示，這個歸一化的預處理工作對最後的結果沒有影響，在後續步驟中也有歸一化的過程，那些過程可以取代這個預處理的歸 一化。詳細請見參考資料[1].

5 梯度的計算(Gradient computation)

簡單地使用一個一維的離散微分模板在一個方向上或者同時在水平和垂直兩個方向上對影象進行處理，更確切地說，這個方法需要使用下面的濾波器核濾除影象中的色彩或變化劇烈的資料作者也嘗試了其他一些更復雜的模板，如3×3 Sobel 模板，或diagonal masks，但是在這個行人檢測的實驗中，這些複雜模板的表現都較差，所以作者的結論是：模板越簡單，效果反而越好。作者也嘗試了在使用微分模板前加入 一個高斯平滑濾波，但是這個高斯平滑濾波的加入使得檢測效果更差，原因是：許多有用的影象資訊是來自變化劇烈的邊緣，而在計算梯度之前加入高斯濾波會把這 些邊緣濾除掉。

6 構建方向的直方圖(creating the orientation histograms)

此步就是為影象的每個細胞單元構建梯度方向直方圖。細胞單元中的每一個畫素點都為某個基於方向的直方圖通道投票。投票是採取加權投票(weighted voting)的方式，即每一票都是帶權值的，這個權值是根據該畫素點的梯度幅度計算出來。可以採用幅值本身或者它的函式來表示這個權值，實際測試表明： 使用幅值來表示權值能獲得最佳的效果，當然，也可以選擇幅值的函式來表示，比如幅值的平方根、幅值的平方、幅值的截斷形式等。細胞單元可以是矩形的，也可以是星形的。直方圖通道是平均分佈在無向(0_{180度)或有向(0}360度)範圍內。
作者發現，採用無向的梯度和9個直方圖通道，能在行人檢測試驗中取得最佳的效果。把細胞單元組合成大的區間由於區域性光照的變化以及前景-背景對比度的變化，使得梯度強度的變化範圍非常大。這就需要對梯度強度做歸一化，作者採取的辦法是：把各個細胞單元組合成大的、空間上連通的區間(blocks)。 這樣以來，HOG描述器就變成了由各區間所有細胞單元的直方圖成分所組成的一個向量。這些區間是互有重疊的，這就意味著：每一個細胞單元的輸出都多次作用 於最終的描述器。區間有兩個主要的幾何形狀:矩形區間(R-HOG)和環形區間(C-HOG)。


R-HOG: R-HOG區間大體上是一些方形的格子，它可以有三個引數來表徵：每個區間中細胞單元的數目、每個細胞單元中畫素點的數目、每個細胞的直方圖通道數目。作者通過實驗表明，行人檢測的最佳引數設定是：3×3細胞 /區間、6×6畫素/細胞、9個直方圖通道。作者還發現，在對直方圖做處理之前，給每個block加一個高斯空域視窗(Gaussian spatial window)是非常必要的，因為這樣可以降低邊緣的周圍畫素點的權重。R- HOG跟SIFT描述器看起來很相似，但他們的不同之處是：R-HOG是在單一尺度下、密集的網格內、沒有對方向排序的情況下被計算出來；而SIFT描述器是在多尺度下、稀疏的影象關鍵點上、對方向排序的情況下被計算出來。補充一點，R-HOG是各區間被組合起來用於對空域資訊進行編碼，而SIFT的各描述器是單獨使用的。
C-HOG: C- HOG區間有兩種不同的形式，它們的區別在於：一個的中心細胞是完整的，一個的中心細胞是被分割的。如上圖所示：作者發現 C-HOG的這兩種形式都能取得相同的效果。C-HOG區間可以用四個引數來表徵：角度盒子的個數、半徑盒子個數、中心盒子的半徑、半徑的伸展因子。通過實驗，對於行人檢測，最佳的引數設定為：4個角度盒子、2個半徑盒子、中心盒子半徑為4個畫素、伸展因子為2。前面提到過，對於R- HOG，中間加一個高斯空域視窗是非常有必要的，但對於C-HOG，這顯得沒有必要。C-HOG看起來很像基於形狀Shape Contexts的方法，但不同之處是：C-HOG的區間中包含的細胞單元有多個orientation channels，而基於形狀上下文的方法僅僅只用到了一個單一的edge presence count。

區間歸一化(Block normalization Schemes)作者採用了四中不同的方法對區間進行歸一化，並對結果進行了比較。引入v表示一個還沒有被歸一 化的向量它包含了給定區間（block）的所有直方圖資訊。| | vk | |表示v的k階範數，這裡的k去1、2。用e表示一個很小的常數。這時，歸一化因子可以表示如下：

L2-norm:L1-norm:L1-sqrt:還 有第四種歸一化方式：L2-Hys，它可以通過先進行L2-norm，對結果進行clipping，然後再重新歸一化得到。作者發現：採用L2- Hys L2-norm 和 L1-sqrt方式所取得的效果是一樣的，L1-norm稍微表現出一點點不可靠性。但是對於沒有被歸一化的資料來說，這四種方法都表現出來顯著的改進。SVM 分類器(SVM classifier)最後一步就是把提取的HOG特徵輸入到SVM分類器中，尋找一個最優超平面作為決策函式。作者採用 的方法是：使用免費的SVMLight軟體包加上HOG分類器來尋找測試影象中的行人。
綜上所述，HoG沒有旋轉和尺度不變性，因此計算量小；而SIFT中每個特徵需要用128維的向量來描述，因此計算量相對很大。而行人檢測中，對於解決Scale-invariant 的問題：將圖片進行不同尺度的縮放，就相當於對模板進行不同尺度scale的縮放.對於解決Rotation-invariant 的問題：建立不同方向的模版（一般取157的）進行匹配.總的來說，就是在不同尺度上的影象進行不同方向的模板(157)匹配，每個點形成一個8方向的梯度描述。SIFT由於其龐大計算量不用與行人檢測，而PCA-SIFT的方法過濾掉很多維度的資訊，只保留20個主分量，因此只適用於行為變化不大的物體檢測。
與其它的特徵描述方法(如：SIFT,　PCA-SIFT,　SURF)相比，HOG有優點是：首先，由於HOG是在影象的區域性方格單元上操作，所以它對影象幾何的和光學的形變都能保持很好的不變性，這兩種形變只會出現在更大的空間領域上。其次，在粗的空域抽樣、精細的方向抽樣以及較強的區域性光學歸一化等條件下，只要行人大體上能夠保持直立的姿勢，可以容許行人有一些細微的肢體動作，這些細微的動作可以被忽略而不影響檢測效果。因此HOG特徵是特別適合於做影象中的人體檢測。

二、原始碼

clc
clear all;
close all;
load('AR_HOG')
%劃分訓練集和測試集
% for i=1:1:120
%    eval(['train_',num2str(i),'=','zeros(20,1824)'';']);
%    %train_set = zeros(20*120,1824);
%    %test_set = zeros(6*120,1824);
%    eval(['test_',num2str(i),'=','zeros(6,1824)'';']); 
% for j=1:1:20
%     eval(['train_',num2str(i),'(:,j)','=','my_lbp(:,26*(i-1)+j)'';']);
% end
% for k=1:1:6
%     eval(['test_',num2str(i),'(:,k)','=','my_lbp(:,26*(i-1)+k+20)'';']);
% end
% end


%% 資料劃分
train_set = zeros(20*120,720);
test_set = zeros(6*120,720);
 train_set = [];
 test_set = [];
for i = 1:120
     train_set = pic_all(:,(i-1)*26+1:(i-1)*26+20);
     test_set = pic_all(:,(i-1)*26+21:(i-1)*26+26);
end
%% 分類
label = zeros(6*120,1);
for i = 1:6*120
    y = test_set(:,i);
    dis = sqrt(sum((train_set - repmat(y,1,20*120)).^2));%計算距離
    local = find(dis == min(dis));
    if mod(local,20)~= 0
        label(i) = fix(local/20)+1;
    else 
        label(i) = local/20;
    end
    
end



%% 計算準確率
a = 1:120;
a = repmat(a,6,1);
real_label = a(:);
dis_label = real_label - label;
acc_rate = length(find(dis_label == 0))/(120*6);








%絕對距離和歐氏距離比較3
% D=0;D1=1;
% for i=1:1:120
%  for j=1:1:19
%   for k=j+1:1:20
% a=eval(['train_',num2str(i),'(:,j)',';']);
% b=eval(['train_',num2str(i),'(:,k)',';'])';
% d=mandist(b,a);
% D=D+d;
% d1=dist(b,a);
% D1=D1+d1;
%   end
%  end
% eval(['ave_d_',num2str(i),'=','D/190'';']);
% eval(['ave_d1_',num2str(i),'=','D1/190'';']);
% D=0;D1=0;
% end
% o=0;f=0;q=0;w=0;
% for i=1:1:120
%  for j=1:1:6
%   for k=1:1:20
%       a=eval(['train_',num2str(i),'(:,k)',';']);
%       b=eval(['test_',num2str(i),'(:,j)',';'])';
%       d=mandist(b,a);
%       d1=dist(b,a);
%   end
%    if d<=1.1*eval(['ave_d_',num2str(i)])
%     f=f+1;
%    else
%        q=q+1;
%    end
%    if d1<=1.1*eval(['ave_d1_',num2str(i)])
%     o=o+1;
%    else
%        w=w+1;
%    end
function  featureVec= hog(img)
img=double(img);
step=8;  %8*8個畫素作為一個cell
[m1 n1]=size(img);
%改變影象尺寸為step的最近整數倍
img=imresize(img,[floor(m1/step)*step,floor(n1/step)*step],'nearest');
[m n]=size(img);
% 2、%伽馬校正
img=sqrt(img);
% 3、求梯度和方向
fy=[-1 0 1];        %定義豎直模板
fx=fy';             %定義水平模板
Iy=imfilter(img,fy,'replicate');    %豎直梯度
Ix=imfilter(img,fx,'replicate');    %水平梯度
Ied=sqrt(Ix.^2+Iy.^2);              %梯度幅值
Iphase=Iy./Ix;              %邊緣斜率，有些為inf,-inf,nan，其中nan需要再處理一下
the=atan(Iphase)*180/3.14159; %求梯度角度


for i=1:m
    for j=1:n
        if(Ix(i,j)>=0&&Iy(i,j)>=0) %第一象限
            the(i,j)=the(i,j);
        elseif(Ix(i,j)<=0&&Iy(i,j)>=0) %第二象限
            the(i,j)=the(i,j)+180;
        elseif(Ix(i,j)<=0&&Iy(i,j)<=0) %第三象限
            the(i,j)=the(i,j)+180;
        elseif(Ix(i,j)>=0&&Iy(i,j)<=0) %第四象限
            the(i,j)=the(i,j)+360;
        end
        if isnan(the(i,j))==1  %0/0會得到nan，如果畫素是nan，重設為0
            the(i,j)=0;
        end
    end
end
the=the+0.000001; %防止角度為

% 4、劃分cell，求cell的直方圖( 1 cell = 8*8 pixel )
%下面是求cell
step=8;                %step*step個畫素作為一個cell
orient=9;               %方向直方圖的方向個數
jiao=360/orient;        %每個方向包含的角度數
Cell=cell(1,1);              %所有的角度直方圖,cell是可以動態增加的，所以先設了一個
ii=1;
jj=1;
for i=1:step:m
    ii=1;
    for j=1:step:n
        Hist1(1:orient)=0;
        for p=1:step
            for q=1:step
                %梯度方向直方圖
                Hist1(ceil(the(i+p-1,j+q-1)/jiao))=Hist1(ceil(the(i+p-1,j+q-1)/jiao))+Ied(i+p-1,j+q-1);
            end
        end
        Cell{ii,jj}=Hist1;       %放入Cell中
        ii=ii+1;
    end
    jj=jj+1;
end
% 5、劃分block，求block的特徵值,使用重疊方式( 1 block = 2*2 cell )
[m n]=size(Cell);
feature=cell(1,(m-1)*(n-1));
for i=1:m-1
    for j=1:n-1
        block=[];
        block=[Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+1,j}(:)' Cell{i+1,j+1}(:)'];
        block=block./sum(block); %歸一化
        feature{(i-1)*(n-1)+j}=block;
    end
end
% 6、影象的HOG特徵值
[m n]=size(feature);
l=2*2*orient;
featureVec=zeros(1,n*l);
for i=1:n
    featureVec((i-1)*l+1:i*l)=feature{i}(:);
end

%  [m n]=size(img);
%  img=sqrt(img);      %伽馬校正
% %下面是求邊緣
% fy=[-1 0 1];        %定義豎直模板
% fx=fy';             %定義水平模板
% Iy=imfilter(img,fy,'replicate');    %豎直邊緣
% Ix=imfilter(img,fx,'replicate');    %水平邊緣
% Ied=sqrt(Ix.^2+Iy.^2);              %邊緣強度
% Iphase=Iy./Ix;              %邊緣斜率，有些為inf,-inf,nan，其中nan需要再處理一下
% 
% 
% %下面是求cell
% step=16;                %step*step個畫素作為一個單元
% orient=9;               %方向直方圖的方向個數
% jiao=360/orient;        %每個方向包含的角度數
% Cell=cell(1,1);              %所有的角度直方圖,cell是可以動態增加的，所以先設了一個
% ii=1;                      
% jj=1;
% for i=1:step:m          %如果處理的m/step不是整數，最好是i=1:step:m-step
%     ii=1;
%     for j=1:step:n      %註釋同上
%          tmpx=Ix(i:i+step-1,j:j+step-1);
%         tmped=Ied(i:i+step-1,j:j+step-1);
%         tmped=tmped/sum(sum(tmped));        %區域性邊緣強度歸一化
%         tmpphase=Iphase(i:i+step-1,j:j+step-1);
%         Hist=zeros(1,orient);               %當前step*step畫素塊統計角度直方圖,就是cell
%         for p=1:step
%             for q=1:step
%                 if isnan(tmpphase(p,q))==1  %0/0會得到nan，如果畫素是nan，重設為0
%                     tmpphase(p,q)=0;
%                 end
%                 ang=atan(tmpphase(p,q));    %atan求的是[-90 90]度之間
%                 ang=mod(ang*180/pi,360);    %全部變正，-90變270
%                 if tmpx(p,q)<0              %根據x方向確定真正的角度
%                     if ang<90               %如果是第一象限
%                         ang=ang+180;        %移到第三象限
%                     end
%                     if ang>270              %如果是第四象限
%                         ang=ang-180;        %移到第二象限
%                     end
%                 end
%                 ang=ang+0.0000001;          %防止ang為0
%                 Hist(ceil(ang/jiao))=Hist(ceil(ang/jiao))+tmped(p,q);   %ceil向上取整，使用邊緣強度加權
%             end
%         end

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423