matlab練習程序（HOG方向梯度直方圖）

HOG（Histogram of Oriented Gradient）方向梯度直方圖，主要用來提取圖像特征，最常用的是結合svm進行行人檢測。

算法流程圖如下（這篇論文上的）：

下面我再結合自己的程序，表述一遍吧：

1.對原圖像gamma校正，img=sqrt(img);

2.求圖像豎直邊緣，水平邊緣，邊緣強度，邊緣斜率。

3.將圖像每16*16（取其他也可以）個像素分到一個cell中。對於256*256的lena來說，就分成了16*16個cell了。

4.對於每個cell求其梯度方向直方圖。通常取9（取其他也可以）個方向（特征），也就是每360/9=40度分到一個方向，方向大小按像素邊緣強度加權。最后歸一化直方圖。

5.每2*2（取其他也可以）個cell合成一個block，所以這里就有（16-1）*（16-1）=225個block。

6.所以每個block中都有2*2*9個特征，一共有225個block，所以總的特征有225*36個。

當然一般HOG特征都不是對整幅圖像取的，而是對圖像中的一個滑動窗口取的。

lena圖：

求得的225*36個特征：

matlab代碼如下：

clear all; close all; clc;

img=double(imread('lena.jpg'));
imshow(img,[]);
[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
        Hist=Hist/sum(Hist);    %方向直方圖歸一化
        Cell{ii,jj}=Hist;       %放入Cell中
        ii=ii+1;                %針對Cell的y坐標循環變量
    end
    jj=jj+1;                    %針對Cell的x坐標循環變量
end

%下面是求feature,2*2個cell合成一個block,沒有顯式的求block
[m n]=size(Cell);
feature=cell(1,(m-1)*(n-1));
for i=1:m-1
   for j=1:n-1           
        f=[];
        f=[f Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+1,j}(:)' Cell{i+1,j+1}(:)'];
        feature{(i-1)*(n-1)+j}=f;
   end
end

%到此結束，feature即為所求
%下面是為了顯示而寫的
l=length(feature);
f=[];
for i=1:l
    f=[f;feature{i}(:)'];  
end 
figure
mesh(f)