【人脸识别】基于HOG特征KNN算法实现人脸识别matlab源码

一、简介

方向梯度直方图(Histogram of Oriented Gradient,HOG)是用于在计算机视觉和图像处理领域，目标检测的特征描述子。该项技术是用来计算图像局部出现的方向梯度次数或信息进行计数。此种方法跟边缘方向直方图、尺度不变特征变换以及形状上下文方法有很多相似。但与它们的不同点是：HOG的计算基于一致空间的密度矩阵来提高准确率。即：在一个网格密集的大小统一的细胞单元上计算，而且为了提高性能，还采用了重叠的局部对比度归一化技术。HoG特征与SVM分类器结合，已经被广泛应用于图像识别中，尤其在行人检测。本节内容介绍HOG相关知识。

1 基本介绍

提出HOG是由Navneet Dalal & Bill Triggs在2005年发表在CVPR中论文[1]。作者给了特征提取的流程图，如下图1所示：

HOG的核心思想是所检测的局部物体外形能够被光强梯度或边缘方向的分布所描述。通过将整幅图像分割成小的连接区域称为cells，每个cell生成一个方向梯度直方图或者cell中pixel的边缘方向，这些直方图的组合可表示出所检测目标的目标）描述子。为改善准确率，局部直方图可以通过计算图像中一个较大区域称为block的光强作为measure被对比标准化，然后用这个measure归一化这个block中的所有cells.这个归一化过程完成了更好的照射/阴影不变性。与其他描述子相比，HOG得到的描述子保持了几何和光学转化不变性除非物体方向改变。而Block与Cells关系如下图所示：

而现在，我们给出作者做的行人检测试验，如下图6所示：

其中，图中(a)表示所有训练图像集的平均梯度；(b)和©分别表示：图像中每一个区间上的最大最大正、负SVM权值；(d)表示一副测试图像；(e)计算完R-HOG后的测试图像；(f)和(g)分别表示被正、负SVM权值加权后的R-HOG图像。

3 算法描述

HOG特征提取方法：首先将图像分成小的连通区域，我们把它叫细胞单元。然后采集细胞单元中各像素点的梯度的或边缘的方向直方图。最后把这些直方图组合起来就可以构成特征描述器。如下图所示。

即：将一个图像：第一步，灰度化即：将图像看做一个x,y,z(灰度)的三维图像；第二步：划分成小cells(2*2)；第三步，计算每个cell中每个pixel的gradient(即orientation)；最后， 统计每个cell的梯度直方图(不同梯度的个数)，即可形成每个cell的descriptor。整个算法具体过程由以下几个部分组成：

4 色彩和伽马归一化(Gamma/Color normalization)

在实际在实际应用中，这一步可以省略，因为作者分别在灰度空间、RGB色彩空间和LAB色彩空间上对图像进行色彩和伽马归一化。而实验结果显示，这个归一化的预处理工作对最后的结果没有影响，在后续步骤中也有归一化的过程，那些过程可以取代这个预处理的归 一化。详细请见参考资料[1].

5 梯度的计算(Gradient computation)

简单地使用一个一维的离散微分模板在一个方向上或者同时在水平和垂直两个方向上对图像进行处理，更确切地说，这个方法需要使用下面的滤波器核滤除图像中的色彩或变化剧烈的数据作者也尝试了其他一些更复杂的模板，如3×3 Sobel 模板，或diagonal masks，但是在这个行人检测的实验中，这些复杂模板的表现都较差，所以作者的结论是：模板越简单，效果反而越好。作者也尝试了在使用微分模板前加入 一个高斯平滑滤波，但是这个高斯平滑滤波的加入使得检测效果更差，原因是：许多有用的图像信息是来自变化剧烈的边缘，而在计算梯度之前加入高斯滤波会把这 些边缘滤除掉。

6 构建方向的直方图(creating the orientation histograms)

此步就是为图像的每个细胞单元构建梯度方向直方图。细胞单元中的每一个像素点都为某个基于方向的直方图通道投票。投票是采取加权投票(weighted voting)的方式，即每一票都是带权值的，这个权值是根据该像素点的梯度幅度计算出来。可以采用幅值本身或者它的函数来表示这个权值，实际测试表明： 使用幅值来表示权值能获得最佳的效果，当然，也可以选择幅值的函数来表示，比如幅值的平方根、幅值的平方、幅值的截断形式等。细胞单元可以是矩形的，也可以是星形的。直方图通道是平均分布在无向(0180度)或有向(0360度)范围内。
作者发现，采用无向的梯度和9个直方图通道，能在行人检测试验中取得最佳的效果。把细胞单元组合成大的区间由于局部光照的变化以及前景-背景对比度的变化，使得梯度强度的变化范围非常大。这就需要对梯度强度做归一化，作者采取的办法是：把各个细胞单元组合成大的、空间上连通的区间(blocks)。 这样以来，HOG描述器就变成了由各区间所有细胞单元的直方图成分所组成的一个向量。这些区间是互有重叠的，这就意味着：每一个细胞单元的输出都多次作用 于最终的描述器。区间有两个主要的几何形状:矩形区间(R-HOG)和环形区间(C-HOG)。


R-HOG: R-HOG区间大体上是一些方形的格子，它可以有三个参数来表征：每个区间中细胞单元的数目、每个细胞单元中像素点的数目、每个细胞的直方图通道数目。作者通过实验表明，行人检测的最佳参数设置是：3×3细胞 /区间、6×6像素/细胞、9个直方图通道。作者还发现，在对直方图做处理之前，给每个block加一个高斯空域窗口(Gaussian spatial window)是非常必要的，因为这样可以降低边缘的周围像素点的权重。R- HOG跟SIFT描述器看起来很相似，但他们的不同之处是：R-HOG是在单一尺度下、密集的网格内、没有对方向排序的情况下被计算出来；而SIFT描述器是在多尺度下、稀疏的图像关键点上、对方向排序的情况下被计算出来。补充一点，R-HOG是各区间被组合起来用于对空域信息进行编码，而SIFT的各描述器是单独使用的。
C-HOG: C- HOG区间有两种不同的形式，它们的区别在于：一个的中心细胞是完整的，一个的中心细胞是被分割的。如上图所示：作者发现 C-HOG的这两种形式都能取得相同的效果。C-HOG区间可以用四个参数来表征：角度盒子的个数、半径盒子个数、中心盒子的半径、半径的伸展因子。通过实验，对于行人检测，最佳的参数设置为：4个角度盒子、2个半径盒子、中心盒子半径为4个像素、伸展因子为2。前面提到过，对于R- HOG，中间加一个高斯空域窗口是非常有必要的，但对于C-HOG，这显得没有必要。C-HOG看起来很像基于形状Shape Contexts的方法，但不同之处是：C-HOG的区间中包含的细胞单元有多个orientation channels，而基于形状上下文的方法仅仅只用到了一个单一的edge presence count。

区间归一化(Block normalization Schemes)作者采用了四中不同的方法对区间进行归一化，并对结果进行了比较。引入v表示一个还没有被归一 化的向量它包含了给定区间（block）的所有直方图信息。| | vk | |表示v的k阶范数，这里的k去1、2。用e表示一个很小的常数。这时，归一化因子可以表示如下：

L2-norm:L1-norm:L1-sqrt:还 有第四种归一化方式：L2-Hys，它可以通过先进行L2-norm，对结果进行clipping，然后再重新归一化得到。作者发现：采用L2- Hys L2-norm 和 L1-sqrt方式所取得的效果是一样的，L1-norm稍微表现出一点点不可靠性。但是对于没有被归一化的数据来说，这四种方法都表现出来显着的改进。SVM 分类器(SVM classifier)最后一步就是把提取的HOG特征输入到SVM分类器中，寻找一个最优超平面作为决策函数。作者采用 的方法是：使用免费的SVMLight软件包加上HOG分类器来寻找测试图像中的行人。
综上所述，HoG没有旋转和尺度不变性，因此计算量小；而SIFT中每个特征需要用128维的向量来描述，因此计算量相对很大。而行人检测中，对于解决Scale-invariant 的问题：将图片进行不同尺度的缩放，就相当于对模板进行不同尺度scale的缩放.对于解决Rotation-invariant 的问题：建立不同方向的模版（一般取157的）进行匹配.总的来说，就是在不同尺度上的图像进行不同方向的模板(157)匹配，每个点形成一个8方向的梯度描述。SIFT由于其庞大计算量不用与行人检测，而PCA-SIFT的方法过滤掉很多维度的信息，只保留20个主分量，因此只适用于行为变化不大的物体检测。
与其它的特征描述方法(如：SIFT,　PCA-SIFT,　SURF)相比，HOG有优点是：首先，由于HOG是在图像的局部方格单元上操作，所以它对图像几何的和光学的形变都能保持很好的不变性，这两种形变只会出现在更大的空间领域上。其次，在粗的空域抽样、精细的方向抽样以及较强的局部光学归一化等条件下，只要行人大体上能够保持直立的姿势，可以容许行人有一些细微的肢体动作，这些细微的动作可以被忽略而不影响检测效果。因此HOG特征是特别适合于做图像中的人体检测。

二、源代码

  

   

    

     

    
    

     

      clc
     
    
   

    

     

    
    

     

      clear all;
     
    
   

    

     

    
    

     

      close all;
     
    
   

    

     

    
    

     

      load('AR_HOG')
     
    
   

    

     

    
    

     

      %
      划分训练集和测试集
     
    
   

    

     

    
    

     

      %
       for i=1:1:120
     
    
   

    

     

    
    

     

      %
       eval(['train_',num2str(i),'=','zeros(20,1824)'';']);
     
    
   

    

     

    
    

     

      %
       %train_set = zeros(20*120,1824);
     
    
   

    

     

    
    

     

      %
       %test_set = zeros(6*120,1824);
     
    
   

    

     

    
    

     

      %
       eval(['test_',num2str(i),'=','zeros(6,1824)'';']); 
     
    
   

    

     

    
    

     

      %
       for j=1:1:20
     
    
   

    

     

    
    

     

      %
       eval(['train_',num2str(i),'(:,j)','=','my_lbp(:,26*(i-1)+j)'';']);
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       for k=1:1:6
     
    
   

    

     

    
    

     

      %
       eval(['test_',num2str(i),'(:,k)','=','my_lbp(:,26*(i-1)+k+20)'';']);
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

       
     
    
   

    

     

    
    

     

      %
      % 数据划分
     
    
   

    

     

    
    

     

      train_set = zeros(20*120,720);
     
    
   

    

     

    
    

     

      test_set = zeros(6*120,720);
     
    
   

    

     

    
    

     

       train_set = [];
     
    
   

    

     

    
    

     

       test_set = [];
     
    
   

    

     

    
    

     

      for i = 1:120
     
    
   

    

     

    
    

     

           train_set = pic_all(:,(i-1)*26+1:(i-1)*26+20);
     
    
   

    

     

    
    

     

           test_set = pic_all(:,(i-1)*26+21:(i-1)*26+26);
     
    
   

    

     

    
    

     

      end
     
    
   

    

     

    
    

     

      %
      % 分类
     
    
   

    

     

    
    

     

      label = zeros(6*120,1);
     
    
   

    

     

    
    

     

      for i = 1:6*120
     
    
   

    

     

    
    

     

          y = test_set(:,i);
     
    
   

    

     

    
    

     

          dis = sqrt(sum((train_set - repmat(y,1,20*120)).^2));%计算距离
     
    
   

    

     

    
    

     

          local = find(dis == min(dis));
     
    
   

    

     

    
    

     

          if mod(local,20)~= 0
     
    
   

    

     

    
    

     

              label(i) = fix(local/20)+1;
     
    
   

    

     

    
    

     

          else 
     
    
   

    

     

    
    

     

              label(i) = local/20;
     
    
   

    

     

    
    

     

          end
     
    
   

    

     

    
    

     
    
     
    
   

    

     

    
    

     

      end
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

       
     
    
   

    

     

    
    

     

      %
      % 计算准确率
     
    
   

    

     

    
    

     

      a = 1:120;
     
    
   

    

     

    
    

     

      a = repmat(a,6,1);
     
    
   

    

     

    
    

     

      real_label = a(:);
     
    
   

    

     

    
    

     

      dis_label = real_label - label;
     
    
   

    

     

    
    

     

      acc_rate = length(find(dis_label == 0))/(120*6);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

       
     
    
   

    

     

    
    

     

      %
      绝对距离和欧氏距离比较3
     
    
   

    

     

    
    

     

      %
       D=0;D1=1;
     
    
   

    

     

    
    

     

      %
       for i=1:1:120
     
    
   

    

     

    
    

     

      %
       for j=1:1:19
     
    
   

    

     

    
    

     

      %
       for k=j+1:1:20
     
    
   

    

     

    
    

     

      %
       a=eval(['train_',num2str(i),'(:,j)',';']);
     
    
   

    

     

    
    

     

      %
       b=eval(['train_',num2str(i),'(:,k)',';'])';
     
    
   

    

     

    
    

     

      %
       d=mandist(b,a);
     
    
   

    

     

    
    

     

      %
       D=D+d;
     
    
   

    

     

    
    

     

      %
       d1=dist(b,a);
     
    
   

    

     

    
    

     

      %
       D1=D1+d1;
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       eval(['ave_d_',num2str(i),'=','D/190'';']);
     
    
   

    

     

    
    

     

      %
       eval(['ave_d1_',num2str(i),'=','D1/190'';']);
     
    
   

    

     

    
    

     

      %
       D=0;D1=0;
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       o=0;f=0;q=0;w=0;
     
    
   

    

     

    
    

     

      %
       for i=1:1:120
     
    
   

    

     

    
    

     

      %
       for j=1:1:6
     
    
   

    

     

    
    

     

      %
       for k=1:1:20
     
    
   

    

     

    
    

     

      %
       a=eval(['train_',num2str(i),'(:,k)',';']);
     
    
   

    

     

    
    

     

      %
       b=eval(['test_',num2str(i),'(:,j)',';'])';
     
    
   

    

     

    
    

     

      %
       d=mandist(b,a);
     
    
   

    

     

    
    

     

      %
       d1=dist(b,a);
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       if d<=1.1*eval(['ave_d_',num2str(i)])
     
    
   

    

     

    
    

     

      %
       f=f+1;
     
    
   

    

     

    
    

     

      %
       else
     
    
   

    

     

    
    

     

      %
       q=q+1;
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       if d1<=1.1*eval(['ave_d1_',num2str(i)])
     
    
   

    

     

    
    

     

      %
       o=o+1;
     
    
   

    

     

    
    

     

      %
       else
     
    
   

    

     

    
    

     

      %
       w=w+1;
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      function  featureVec= hog(img)
     
    
   

    

     

    
    

     

      img=double(img);
     
    
   

    

     

    
    

     

      step=8;  %8*8个像素作为一个cell
     
    
   

    

     

    
    

     

      [m1 n1]=size(img);
     
    
   

    

     

    
    

     

      %
      改变图像尺寸为step的最近整数倍
     
    
   

    

     

    
    

     

      img=imresize(img,[floor(m1/step)*step,floor(n1/step)*step],'nearest');
     
    
   

    

     

    
    

     

      [m n]=size(img);
     
    
   

    

     

    
    

     

      %
       2、%伽马校正
     
    
   

    

     

    
    

     

      img=sqrt(img);
     
    
   

    

     

    
    

     

      %
       3、求梯度和方向
     
    
   

    

     

    
    

     

      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需要再处理一下
     
    
   

    

     

    
    

     

      the=atan(Iphase)*180/3.14159; %求梯度角度
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

      for i=1:m
     
    
   

    

     

    
    

     

          for j=1:n
     
    
   

    

     

    
    

     

              if(Ix(i,j)>=0&&Iy(i,j)>=0) %第一象限
     
    
   

    

     

    
    

     

                  the(i,j)=the(i,j);
     
    
   

    

     

    
    

     

              elseif(Ix(i,j)<=0&&Iy(i,j)>=0) %第二象限
     
    
   

    

     

    
    

     

                  the(i,j)=the(i,j)+180;
     
    
   

    

     

    
    

     

              elseif(Ix(i,j)<=0&&Iy(i,j)<=0) %第三象限
     
    
   

    

     

    
    

     

                  the(i,j)=the(i,j)+180;
     
    
   

    

     

    
    

     

              elseif(Ix(i,j)>=0&&Iy(i,j)<=0) %第四象限
     
    
   

    

     

    
    

     

                  the(i,j)=the(i,j)+360;
     
    
   

    

     

    
    

     

              end
     
    
   

    

     

    
    

     

              if isnan(the(i,j))==1  %0/0会得到nan，如果像素是nan，重设为0
     
    
   

    

     

    
    

     

                  the(i,j)=0;
     
    
   

    

     

    
    

     

              end
     
    
   

    

     

    
    

     

          end
     
    
   

    

     

    
    

     

      end
     
    
   

    

     

    
    

     

      the=the+0.000001; %防止角度为
     
    
   

    

     

    
    

     

      %
       4、划分cell，求cell的直方图( 1 cell = 8*8 pixel )
     
    
   

    

     

    
    

     

      %
      下面是求cell
     
    
   

    

     

    
    

     

      step=8;                %step*step个像素作为一个cell
     
    
   

    

     

    
    

     

      orient=9;               %方向直方图的方向个数
     
    
   

    

     

    
    

     

      jiao=360/orient;        %每个方向包含的角度数
     
    
   

    

     

    
    

     

      Cell=cell(1,1);              %所有的角度直方图,cell是可以动态增加的，所以先设了一个
     
    
   

    

     

    
    

     

      ii=1;
     
    
   

    

     

    
    

     

      jj=1;
     
    
   

    

     

    
    

     

      for i=1:step:m
     
    
   

    

     

    
    

     

          ii=1;
     
    
   

    

     

    
    

     

          for j=1:step:n
     
    
   

    

     

    
    

     

              Hist1(1:orient)=0;
     
    
   

    

     

    
    

     

              for p=1:step
     
    
   

    

     

    
    

     

                  for q=1:step
     
    
   

    

     

    
    

     

                      %梯度方向直方图
     
    
   

    

     

    
    

     

                      Hist1(ceil(the(i+p-1,j+q-1)/jiao))=Hist1(ceil(the(i+p-1,j+q-1)/jiao))+Ied(i+p-1,j+q-1);
     
    
   

    

     

    
    

     

                  end
     
    
   

    

     

    
    

     

              end
     
    
   

    

     

    
    

     

              Cell{ii,jj}=Hist1;       %放入Cell中
     
    
   

    

     

    
    

     

              ii=ii+1;
     
    
   

    

     

    
    

     

          end
     
    
   

    

     

    
    

     

          jj=jj+1;
     
    
   

    

     

    
    

     

      end
     
    
   

    

     

    
    

     

      %
       5、划分block，求block的特征值,使用重叠方式( 1 block = 2*2 cell )
     
    
   

    

     

    
    

     

      [m n]=size(Cell);
     
    
   

    

     

    
    

     

      feature=cell(1,(m-1)*(n-1));
     
    
   

    

     

    
    

     

      for i=1:m-1
     
    
   

    

     

    
    

     

          for j=1:n-1
     
    
   

    

     

    
    

     

              block=[];
     
    
   

    

     

    
    

     

              block=[Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+1,j}(:)' Cell{i+1,j+1}(:)'];
     
    
   

    

     

    
    

     

              block=block./sum(block); %归一化
     
    
   

    

     

    
    

     

              feature{(i-1)*(n-1)+j}=block;
     
    
   

    

     

    
    

     

          end
     
    
   

    

     

    
    

     

      end
     
    
   

    

     

    
    

     

      %
       6、图像的HOG特征值
     
    
   

    

     

    
    

     

      [m n]=size(feature);
     
    
   

    

     

    
    

     

      l=2*2*orient;
     
    
   

    

     

    
    

     

      featureVec=zeros(1,n*l);
     
    
   

    

     

    
    

     

      for i=1:n
     
    
   

    

     

    
    

     

          featureVec((i-1)*l+1:i*l)=feature{i}(:);
     
    
   

    

     

    
    

     

      end
     
    
   

    

     

    
    

     

      %
       [m n]=size(img);
     
    
   

    

     

    
    

     

      %
       img=sqrt(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需要再处理一下
     
    
   

    

     

    
    

     

      %
       
     
    
   

    

     

    
    

     

      %
       
     
    
   

    

     

    
    

     

      %
       %下面是求cell
     
    
   

    

     

    
    

     

      %
       step=16; %step*step个像素作为一个单元
     
    
   

    

     

    
    

     

      %
       orient=9; %方向直方图的方向个数
     
    
   

    

     

    
    

     

      %
       jiao=360/orient; %每个方向包含的角度数
     
    
   

    

     

    
    

     

      %
       Cell=cell(1,1); %所有的角度直方图,cell是可以动态增加的，所以先设了一个
     
    
   

    

     

    
    

     

      %
       ii=1; 
     
    
   

    

     

    
    

     

      %
       jj=1;
     
    
   

    

     

    
    

     

      %
       for i=1:step:m %如果处理的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
     
    
   

    

     

    
    

     

      %
       for p=1:step
     
    
   

    

     

    
    

     

      %
       for q=1:step
     
    
   

    

     

    
    

     

      %
       if isnan(tmpphase(p,q))==1 %0/0会得到nan，如果像素是nan，重设为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/jiao))=Hist(ceil(ang/jiao))+tmped(p,q); %ceil向上取整，使用边缘强度加权
     
    
   

    

     

    
    

     

      %
       end
     
    
   

    

     

    
    

     

      %
       end
     
    
  

三、运行结果

四、备注

完整代码添加QQ1575304183