clear all; close all; clc;
img=imread('3.jpg');
figure,subplot(2,3,1);imshow(img);title('原始彩色图像')
img1= rgb2gray(img);%转化为灰度图像
subplot(2,3,2);imshow(img1);title('灰度图像')
[m n]=size(img1);
img2=double(img1);
img3=sqrt(img2); %伽马校正,img2为一个矩阵，sqrt(img2)可以对矩阵的每一个元素进行开方操作
% subplot(2,3,3);imshow(img3);title('伽马校正后图像')
%下面是求边缘,计算出图像处每一点处的梯度值
fy=[-1 0 1]; %定义竖直模板
fx=fy'; %定义水平模板
Iy=imfilter(img3,fy,'replicate'); %竖直边缘
%可用以下代码实现梯度值得计算，以代替imfilter函数
Ix=imfilter(img3,fx,'replicate'); %水平边缘
Ied=sqrt(Ix.^2+Iy.^2); %边缘强度
Iphase=Iy./Ix; %求出某一点的斜率，y/x,为这个点的正切值，边缘斜率，有些为inf,-inf,nan，其中nan需要再处理一下
%下面是求cell
step=8; %step*step个像素作为一个单元，cell一般是8*8的，步长一般选为8个像素
orient=9; %方向直方图的方向个数
jiao=360/orient; %每个方向包含的角度数,jiao=40
Cell=cell(1,1); %所有的角度直方图,cell是可以动态增加的，所以先设了一个1行1列的空矩阵
ii=1;
jj=1;
for i=1:step:m-step %如果处理的m/step不是整数，最好是i=1:step:m-step，若刚好是step的整数倍，就选择到m，这个跟要检测的图片尺寸有关系
ii=1;
for j=1:step:n-step %注释同上
tmpx=Ix(i:i+step-1,j:j+step-1); % Ix记录了每一点的x方向梯度,在Ix中取出一个step*step，即6*6大小的块出来，
tmped=Ied(i:i+step-1,j:j+step-1);%Ied记录了每一点的梯度方向（含x和y),在Ied中取出一个step*step，即6*6大小的块出来
tmped1=tmped/sum(sum(tmped)); %局部边缘强度归一化,将取出来的6*6块中的每一个元素值除以这一块的总和，对一个矩阵来说，sum(a)等于
%是对a中各列取和，得到一个行向量,若a是一个行向量，sum(a)等于是将a中各个元素相加
%例如a=a=[ 1 3 7;8 6 9;7 3 4],则sum(a) = 16 12 20,sum(sum(a))=48,
tmpphase=Iphase(i:i+step-1,j:j+step-1);
Hist=zeros(1,orient); %Hist为1行9列的值全为0的行向量。当前step*step像素块统计角度直方图,就是cell
for p=1:step
for q=1:step
if isnan(tmpphase(p,q))==1 %0/0会得到nan，如果像素是nan，重设为0，isnan(a)函数就是判断a是否是一个有限大小的数，
%NaN的意思就是不是一个数，我们可以理解为无穷数。0/0就是NAN，
% NaN或者nan都是“非数”的意思，“0/0”、“∞/∞”、“0*∞”都会产生这种结果
tmpphase(p,q)=0;
end
ang1=atan(tmpphase(p,q)); %atan求的是弧度值，
ang2=mod(ang1*180/pi,360); %弧度转换成角度，MOD(number,divisor)函数是用于返回两数相除的余数，返回结果的符号与除数（divisor）的符号相同。
%取值规律,若两个数符号不同，先将两个整数看作是正数，再作除法运算，不能整除时，其值=除数×(整商+1)-被除数，比如
%mod(-45,360)=315,mod(-1,360)=359,若ang本身为正数，即第一象限角，则取模(mod)后，还是第一象限角，若ang为负值，
%其为第三象限的角，与360取模运算后，变成第4象限的角。
%也就是此项运算目的将所有角度换算成正角度
if tmpx(p,q)<0 %如果所判断点的x方向梯度为负，即顺着x方向像素值减小
if ang2<90 %如果是第一象限
ang2=ang2+180; %移到第三象限
end
if ang2>270 %如果是第四象限
ang2=ang2-180; %移到第二象限
end
end
ang2=ang2+0.0000001; %防止ang为0
Hist(ceil(ang2/jiao))=Hist(ceil(ang2/jiao))+tmped1(p,q); %ceil函数为向上取整，即ceil(3.1)=4,使用边缘强度加权，tmped(p,q)是一个6*6的小cell,里面是每一个像素的梯度值
%jiao=40,先判断cell里面即6*6像素块中每个点根据它自身的正切值判断在那个区域（9个区域），
%再用这点自身的梯度
