轉摘網址：http://www.cnblogs.com/tiandsp/archive/2013/05/24/3097503.html
Hog參考網址：http://www.cnblogs.com/tornadomeet/archive/2012/08/15/2640754.html
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個特征：
1 function Hog 2 clear all; close all; clc; 3 4 img=double(imread('lena.jpg')); 5 imshow(img,[]); 6 [m n]=size(img); 7 8 img=sqrt(img); %伽馬校正 9 10 %下面是求邊緣 11 fy=[-1 0 1]; %定義豎直模板 12 fx=fy'; %定義水平模板 13 Iy=imfilter(img,fy,'replicate'); %豎直邊緣 14 Ix=imfilter(img,fx,'replicate'); %水平邊緣 15 Ied=sqrt(Ix.^2+Iy.^2); %邊緣強度 16 Iphase=Iy./Ix; %邊緣斜率，有些為inf,-inf,nan，其中nan需要再處理一下 17 18 19 %下面是求cell 20 step=16; %step*step個像素作為一個單元 21 orient=9; %方向直方圖的方向個數 22 jiao=360/orient; %每個方向包含的角度數 23 Cell=cell(1,1); %所有的角度直方圖,cell是可以動態增加的，所以先設了一個 24 ii=1; 25 jj=1; 26 for i=1:step:m %如果處理的m/step不是整數，最好是i=1:step:m-step 27 ii=1; 28 for j=1:step:n %注釋同上 29 tmpx=Ix(i:i+step-1,j:j+step-1); 30 tmped=Ied(i:i+step-1,j:j+step-1); 31 tmped=tmped/sum(sum(tmped)); %局部邊緣強度歸一化 32 tmpphase=Iphase(i:i+step-1,j:j+step-1); 33 Hist=zeros(1,orient); %當前step*step像素塊統計角度直方圖,就是cell 34 for p=1:step 35 for q=1:step 36 if isnan(tmpphase(p,q))==1 %0/0會得到nan，如果像素是nan，重設為0 37 tmpphase(p,q)=0; 38 end 39 ang=atan(tmpphase(p,q)); %atan求的是[-90 90]度之間 40 ang=mod(ang*180/pi,360); %全部變正，-90變270 41 if tmpx(p,q)<0 %根據x方向確定真正的角度 42 if ang<90 %如果是第一象限 43 ang=ang+180; %移到第三象限 44 end 45 if ang>270 %如果是第四象限 46 ang=ang-180; %移到第二象限 47 end 48 end 49 ang=ang+0.0000001; %防止ang為0 50 Hist(ceil(ang/jiao))=Hist(ceil(ang/jiao))+tmped(p,q); %ceil向上取整，使用邊緣強度加權 51 end 52 end 53 Hist=Hist/sum(Hist); %方向直方圖歸一化 54 Cell{ii,jj}=Hist; %放入Cell中 55 ii=ii+1; %針對Cell的y坐標循環變量 56 end 57 jj=jj+1; %針對Cell的x坐標循環變量 58 end 59 60 %下面是求feature,2*2個cell合成一個block,沒有顯式的求block 61 [m n]=size(Cell); 62 feature=cell(1,(m-1)*(n-1)); 63 for i=1:m-1 64 for j=1:n-1 65 f=[]; 66 f=[f Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+1,j}(:)' Cell{i+1,j+1}(:)']; 67 feature{(i-1)*(n-1)+j}=f; 68 end 69 end 70 71 %到此結束，feature即為所求 72 %下面是為了顯示而寫的 73 l=length(feature); 74 f=[]; 75 for i=1:l 76 f=[f;feature{i}(:)']; 77 end 78 figure 79 mesh(f)