HOG算法基础

实现思路步骤：

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:matlena.jpg');%图片位置
%获取图像，尺寸，并将图像resize成step的最近整数倍
img=double(img);
figure;imshow(img,[]);%显示图像
step=8;      %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=[-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需要再处理一下
figure; imshow(Ied,[]);             %显示梯度提取后的值
%求cell
orient=9;                  %方向直方图的方向个数
angular=360/orient;        %每个方向包含的角度数，划分角度区间，0到40度一个区间...
Cell=cell(1,1);            %所有的角度直方图,cell是可以动态增加的，先设了一个
%% 开始获取orient个方向的特征向量
ii=1;                      
jj=1;
for i=1:step:m-step         %如果处理的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

        %% 统计一个cell里面的梯度信息
         for p=1:step
            for q=1:step
                %%  判断是不是一个数字True for Not-a-Number.如果不是一个数字，就归零
                if isnan(tmpphase(p,q))==1  %因为会遇到0/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/angular)) = Hist(ceil(ang/angular))+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
[m2, n2]=size(Cell);
feature=cell(1,(m2-1)*(n2-1));
for i=1:m2-1
   for j=1:n2-1           
        f=[];
        f=[f Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+1,j}(:)' Cell{i+1,j+1}(:)'];
        feature{(i-1)*(n2-1)+j}=f;
   end
end
%以上就可以得到我们需要的所有的HOG特征了，以下部分主要是为了显示
l=length(feature);
f=[];
for i=1: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:UsersdingDesktopTest.jpg');
img=rgb2gray(img); %简单起见，彩图转灰度图。后续可以改进。
img=imresize(img, [128 64]);
img=double(img);

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

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

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

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

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

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

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


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