clear?all;?close?all;?clc;

img=double（imread（‘C:\\Users\\lihaibiao\\Desktop\\test\\David0.jpg‘））;
imshow（img[]）;
[m?n]=size（img）;

img=sqrt（img）;??????%伽馬校正

%下面是求邊緣
fy=[-1?0?1];????????%定義豎直模板
fx=fy‘;?????????????%定義水平模板
Iy=imfilter（imgfy‘replicate‘）;????%豎直邊緣
Ix=imfilter（imgfx‘replicate‘）;????%水平邊緣
Ied=sqrt（Ix.^2+Iy.^2）;??????????????%邊緣強(qiáng)度
Iphase=Iy./Ix;??????????????%邊緣斜率，有些為inf-infnan，其中nan需要再處理一下


%下面是求cell
step=16;????????????????%step*step個像素作為一個單元
orient=9;???????????????%方向直方圖的方向個數(shù)
jiao=360/orient;????????%每個方向包含的角度數(shù)
Cell=cell（11）;??????????????%所有的角度直方圖cell是可以動態(tài)增加的，所以先設(shè)了一個
ii=1;??????????????????????
jj=1;
for?i=1:step:m??????????%如果處理的m/step不是整數(shù)，最好是i=1:step:m-step
????ii=1;
????for?j=1:step:n??????%注釋同上
????????tmpx=Ix（i:i+step-1j:j+step-1）;
????????tmped=Ied（i:i+step-1j:j+step-1）;
????????tmped=tmped/sum（sum（tmped））;????????%局部邊緣強(qiáng)度歸一化
????????tmpphase=Iphase（i:i+step-1j:j+step-1）;
????????Hist=zeros（1orient）;???????????????%當(dāng)前step*step像素塊統(tǒng)計角度直方圖就是cell
????????for?p=1:step
????????????for?q=1:step
????????????????if?isnan（tmpphase（pq））==1??%0/0會得到nan，如果像素是nan，重設(shè)為0
????????????????????tmpphase（pq）=0;
????????????????end
????????????????ang