实现思路步骤：

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特征都不是对整幅图像取的，而是对图像中的一个滑动窗口取的。

形象化的用一个流程图显示：

matlab实现代码：参考别人的修改的

clear;clc;
img=imread('E:\mat\lena.jpg');%图片位置
%获取图像，尺寸，并将图像resize成step的最近整数倍
img=double(img);
figure;imshow(img,[]);%显示图像
step=; %step*step个像素作为一个cell
[m1 ,n1]=size(img);%获取图像尺寸
img=imresize(img,[floor(m1/step)*step,floor(n1/step)*step],'nearest');
[m,n]=size(img);
%伽马校正
img=sqrt(img);
figure; imshow(img,[]);%显示图像
%滤波，求梯度
fy=[- ]; %定义竖直模板
fx=fy'; %定义水平模板
Iy=imfilter(img,fy,'replicate'); %竖直边缘
Ix=imfilter(img,fx,'replicate'); %水平边缘
Ied=sqrt(Ix.^+Iy.^); %边缘强度 求梯度的长度
Iphase=Iy./Ix; %边缘斜率，有些为inf,-inf,nan，其中nan需要再处理一下
figure; imshow(Ied,[]); %显示梯度提取后的值
%求cell
orient=; %方向直方图的方向个数
angular=/orient; %每个方向包含的角度数，划分角度区间，0到40度一个区间…
Cell=cell(,); %所有的角度直方图,cell是可以动态增加的，先设了一个
%% 开始获取orient个方向的特征向量
ii=;
jj=;
for i=:step:m-step %如果处理的m/step不是整数，最好是i=:step:m-step
ii=;
for j=:step:n
tmpx=Ix(i:i+step-,j:j+step-); %水平
tmped=Ied(i:i+step-,j:j+step-); %边缘强度
tmped=tmped/sum(sum(tmped)); %% 局部边缘强度归一化
tmpphase=Iphase(i:i+step-,j:j+step-);%% 边缘斜率局部提取
Hist=zeros(,orient); %% 创建直方图
%当前step*step像素块统计角度直方图,就是cell

    %% 统计一个cell里面的梯度信息  
     for p=:step  
        for q=:step  
            %%  判断是不是一个数字True for Not-a-Number.如果不是一个数字，就归零  
            if isnan(tmpphase(p,q))==  %因为会遇到0/0的情况  
                tmpphase(p,q)=;  
            end  
            %% 进行区间的划分  
            ang=atan(tmpphase(p,q));    %atan求的是\[- \]度之间  
            ang=mod(ang\*/pi,);    %全部变正，-90变270  
            if tmpx(p,q)<              %根据x方向确定真正的角度  
                if ang<               %如果是第一象限  
                    ang=ang+;        %移到第三象限  
                end  
                if ang>              %如果是第四象限  
                    ang=ang-;        %移到第二象限  
                end  
            end  
            ang=ang+0.0000001;          %防止ang为0  
            Hist(ceil(ang/angular)) = Hist(ceil(ang/angular))+tmped(p,q);   %ceil向上取整，使用边缘强度加权  
        end  
     end  
    %% 方向直方图归一化  
    Hist=Hist/sum(Hist);  
    Cell{ii,jj}=Hist;       %放入Cell中  
    ii=ii+;                %针对Cell的y坐标循环变量  
end  
jj=jj+;                    %针对Cell的x坐标循环变量  

end

%% 下面是求feature,*2个cell合成一个block,没有显式的求block
[m2, n2]=size(Cell);
feature=cell(,(m2-)*(n2-));
for i=:m2-
for j=:n2-
f=[];
f=[f Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+,j}(:)' Cell{i+1,j+1}(:)'];
feature{(i-)*(n2-)+j}=f;
end
end
%以上就可以得到我们需要的所有的HOG特征了，以下部分主要是为了显示
l=length(feature);
f=[];
for i=:l
f=[f;feature{i}(:)'];
end
figure
mesh(f)

当然看文献发现我们还可以对得到的HOG特征进行PCA降维处理来提高计算速率的同时，降低噪声的影响。----有点矛盾吧，HOG我感觉实际上是采用不断划分cell和block的过程增加描述图像的特征，而PCA则是降低图像特征描述的维数。

附一个基础整理的mind图：

亲测可用的程序

clear all; close all; clc;

%img=double(imread('lena.jpg'));
%img=imread('man.png');
img=imread('C:\Users\ding\Desktop\Test.jpg');
img=rgb2gray(img); %简单起见，彩图转灰度图。后续可以改进。
img=imresize(img, [ ]);
img=double(img);

[h, w, ~] = size(img);

%下面是求cell
cell_size=; %step*step个像素作为一个cell. cell_size=pixels_per_cell
orient=; %方向直方图包含的方向数
angle_range=/orient; %每个方向包含的角度数

h=round(h/cell_size)*cell_size;
w=round(w/cell_size)*cell_size;
img=img(:h,:w,:);

img = sqrt(img); %伽马校正。J=AI^r 此处取A=,r=0.5

% 下面是求边缘
fy=[- ]; %定义竖直模版
fx=fy'; %定义水平模版

Gy=imfilter(img, fy, 'replicate'); %竖直梯度
Gx=imfilter(img, fx, 'replicate'); %水平梯度
Gmag=sqrt(Gx.^+Gy.^); %梯度幅值

%为每个cell计算其decriptor(梯度方向直方图，即一个1*orient规格的向量）
cell_descriptors=zeros(orient, h/cell_size, w/cell_size);
idx_y=;
for y=:cell_size:h
idx_x=;
for x=:cell_size:w
tmpx=Gx(y:y+cell_size-, x:x+cell_size-);
tmpy=Gy(y:y+cell_size-, x:x+cell_size-);
tmped=Gmag(y:y+cell_size-,x:x+cell_size-);
tmped=tmped/sum(sum(tmped)); %每个cell的局部边缘强度归一化
cell_hist=zeros(, orient); %当前cell_size*cell_size像素统计角度直方图，就是cell
for p=:cell_size
for q=:cell_size
ang=atan2(tmpy(p,q), tmpx(p,q)); %atan2返回的是[-pi,pi]之间的弧度值
ang=mod(ang*/pi, ); %先转角度，再划归到[,)之间。因为mod的参数现在不是整数，因此会大于179.
ang=ang+0.0000001; %防止ang为0

            bin\_id = ceil(ang/angle\_range);%得到的bin\_id \\in \[,\]  
            cell\_hist(bin\_id)=cell\_hist(bin\_id)+tmped(p,q); %ceil向上取整，使用边缘强度加权。此处根据梯度方向进行vote，权值为梯度幅值  
        end  
    end  
    cell\_descriptors(:,idx\_y,idx\_x) = cell\_hist;  
    idx\_x = idx\_x + ;  
end  
idx\_y = idx\_y + ;  

end

%下面是计算feature,block_size*block_size个cell合成一个block
%比如block_size取2
[sw, h, w]=size(cell_descriptors);
block_size=; %cells_per_block=，即每个block_size=*=16像素
%d得到横竖方向数分别有几个block，用于计算总的特征数量
stride=;
h_max=floor((h-block_size)/stride)+;
w_max=floor((w-block_size)/stride)+;
block_descriptors=zeros(block_size*block_size*orient, h_max, w_max);
for i=:h_max
for j=:w_max
blk_mat=cell_descriptors(:,i:i+block_size-, j:j+block_size-);%得到每个block的特征值
%%L2范数归一化就是向量中每个元素除以向量的L2范数
for k=:
A=reshape(blk_mat(k,:,:),,);
normed_blk_mat(k)=norm(A,);
% blk_mat1=
end
normed_blk_mat1=norm(normed_blk_mat,);
normed_blk_mat2=blk_mat/normed_blk_mat1;

    reshaped\_blk\_mat=reshape(normed\_blk\_mat2, \[ block\_size\*block\_size\*orient\]);  
    block\_descriptors(:,i,j)=reshaped\_blk\_mat;  
end  

end

%将block_descriptors进行拼接，得到final_descriptor
[d1,d2,d3]=size(block_descriptors);
dimensions=d1*d2*d3;
final_descriptor=zeros(, dimensions);
k=;
for i=:d2
for j=:d3
final_descriptor(k:k+d1-)=block_descriptors(:,i,j);
k=k+d1;
end
end