搜索
clearall;closeall;clc;img=imread('3.jpg');figure,subplot(2,3,1);imshow(img);title('原始彩色图像')img1=rgb2gray(img);%转化为灰度图像subplot(2,3,2);imshow(img1);title('灰度图像')[mn]=size(img1);img2=double(img1);img3=sqrt(img2);%伽马校正,img2为一个矩阵，sqrt(img2)可以对矩阵的每一个元素进行开方操作%subplot(2,3,3);imshow(img3);title('伽马校正后图像')%下面是求边缘,计算出图像处每一点处的梯度值fy=[-101];%定义竖直模板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需要再处理一下%下面是求cellstep=8;%step*step个像素作为一个单元，cell一般是8*8的，步长一般选为8个像素orient=9;%方向直方图的方向个数jiao=360/orient;%每个方向包含的角度数,jiao=40Cell=cell(1,1);%所有的角度直方图,cell是可以动态增加的，所以先设了一个1行1列的空矩阵ii=1;jj=1;fori=1:step:m-step%如果处理的m/step不是整数，最好是i=1:step:m-step，若刚好是step的整数倍，就选择到m，这个跟要检测的图片尺寸有关系ii=1;forj=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=[137;869;734],则sum(a)=161220,sum(sum(a))=48,tmpphase=Iphase(i:i+step-1,j:j+step-1);Hist=zeros(1,orient);%Hist为1行9列的值全为0的行向量。当前step*step像素块统计角度直方图,就是cellforp=1:stepforq=1:stepifisnan(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;endang1=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象限的角。%也就是此项运算目的将所有角度换算成正角度iftmpx(p,q)0%如果所判断点的x方向梯度为负，即顺着x方向像素值减小ifang290%如果是第一象限ang2=ang2+180;%移到第三象限endifang2270%如果是第四象限ang2=ang2-180;%移到第二象限endendang2=ang2+0.0000001;%防止ang为0Hist(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个区域），%再用这点自身的梯度（Ied=sqrt(Ix.^2+Iy.^2);tmped=Ied(i:i+step-1,j:j+step-1)，tmped=tmped/sum(sum(tmped))值加上去加权一下，%hist的9个值里面，每个值是,endend%这个循环结束后，会将8*8的cell中所有点8*8=64点处的角度算完，并判断每个点会在9个区的哪个区，在那个区就放进hist中的某个值中，%Hist=Hist/sum(Hist);%方向直方图归一化，这一步可以没有，因为是组成block以后再进行归一化就可以Cell{ii,jj}=Hist;%一个8*8的cell中所有点8*8=64个点处的角度算完后，放入Cell中,将每个点的hist值记录下来放进cell里面ii=ii+1;%这个循环中，原图像的j=1:step:n-step,即让j进行循环，即保持横坐标不动，让纵坐标在垂直方向遍历，若垂直方向有128个像素，%步长为8的话，就会遍历15个窗口，即ii最终值为15，cell(1,1)...cell(15,1)共15个cell,cell中的每一个元素记录了一个cell的方向特征，垂直方向循环完后，这样针对Cell的y坐标循环变量endjj=jj+1;%现在开始在水平方向移动step个步长，比如第一步后jj=2,i=8,这是会又会产生15个cell,记为cell(1,2),cell(2,2)...cell(15,2)%后面继续，就会得到所有的cell，最后一个应该是cell(15,m-step），这样一共有（128/8-1）*（64/8-1）=15*7=105个cell,针对Cell的x坐标循环变量,最后生成一个82*54的cellend%下面是求feature,2*2个cell合成一个block,没有显式的求block[m1n1]=size(Cell);feature=cell(1,(m1-1)*(n1-1));fori=1:m1-1forj=1:n1-1f=[];f=[fCell{i,j}(:)'Cell{i,j+1}(:)'Cell{i+1,j}(:)'Cell{i+1,j+1}(:)'];f=f./sum(f);%归一化feature{(i-1)*(n1-1)+j}=f;endend%到此结束，feature即为所求,即一个图片中的所有hog计算完毕%下面是为了显示而写的l=length(feature);f=[];fori=1:lf=[f;feature{i}(:)'];endsubplot(2,3,3);mesh(f)