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)

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