经过Canny算子边界提取后的图像，对于边界点求方向角，公式如下：
其中g(i,j)为经过高斯滤波后图像的灰度。以5度为单位，将360度的角度空间量化为72级，计算每个边界点处法向量的方向角分别落在这72级空间中的频率，这样便得到了边界方向直方图H。 边界方向不受图像中对象的位置的影响，为了达到不受图像缩放的影响，对得到的边界方向直方图进行归一化处理：H=H/S。其中H[I]为边界方向直方图，S为图像的面积。
边界方向直方图很显然受到图像旋转的影响，为了减少图像旋转对于直方图的影响，我们对直方图进行平滑：
参数k决定平滑度。

源码在下面：

% 作者：丕子 472308749@qq.com
% sigma可取 2.5
function eo=edge_o(p,sigma)

eo='';

 a = imread(p);
 if isrgb(a)%判断是否是RGB
 a=rgb2gray(a);
 end
 a= double(a);  % 将图像像素数据转换为double类型

[m,n] = size(a);
e = repmat(logical(uint8(0)),m,n);  % 生成初始矩阵

OffGate = 0.0001;
Perc = 0.7;
Th = 0.4;

pw = 1:30;             % hardcoding. But it's large enough if sigma isn't too large
sig2 = sigma * sigma;  % 方差

width = max(find(exp(-(pw.*pw)/(2*sig2)) > OffGate));  % 寻找截断点

t = (-width:width);
len = 2*width+1;
t3=[t-0.5;t;t+0.5];

dgau = (-t.*exp(-(t.*t)/(2*sig2))/sig2).';              % 一阶高斯函数的导数

ra = size(a,1);  % 图像行数
ca = size(a,2);  % 图像列数
ay = 255*a;
ax = 255*a';

ax = conv2(ax,dgau,'same').';   % 高斯平滑滤波后的图像的x方向梯度
ay = conv2(ay,dgau,'same');     % 高斯平滑滤波后的图像的y方向梯度
mag = sqrt((ax.*ax)+(ay.*ay));  % 每个像素点的梯度强度

magmax = max(mag(:));
if magmax>0
 mag = mag/magmax;  % 归一化
end

[counts,x] = imhist(mag,64);    % 直方图统计
high = min(find(cumsum(counts)>Perc*m*n))/64;
low = Th*high;

thresh = [low high];   % 根据直方图统计确定上下限

% Nonmax-suppression
idxStrong = [];
for dir = 1:4
 Localmax = Findlocalmax(dir,ax,ay,mag);
 idxWeak = Localmax(mag(Localmax)>low);
 e(idxWeak) = 1;%idxWeak就是边界点
 idxStrong = [idxStrong; idxWeak(mag(idxWeak)>high)];
end

rstrong = rem(idxStrong-1,m)+1;
cstrong = floor((idxStrong-1)/m)+1;
e = bwselect(e,cstrong,rstrong,8);  %
e = bwmorph(e,'thin',1);   % 使用形态学的方法把边缘变细

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算边界方向-lipiji
H=zeros(1,72);
ii=1;
HA=zeros(1,72);
for jj=-180:5:180
 HA(ii)=jj;
 ii=ii+1;
end

count1=0;
for i=1:ra*ca
 if e(i)==1
 theta =atan2(ay(i),ax(i));
 theta=round(theta*180/3.141592653);%变成角度制

 mina=1000;
 for ii=1:72
 if abs(theta-HA(ii))<mina
 mina=abs(theta-HA(ii));
 mini=ii;
 end
 end
 H(mini)=H(mini)+1;
 count1=count1+1;
 end

end
 H=H/count1;%归一化，保证尺度不变性，count1为边界总长度
 k=2;%平滑系数。决定平滑度

 for i=(1+k):(72-k)
 summ=0;
 for j=(i-k):(i+k)
 summ=summ+H(j);
 end
 H(i)=summ/(2*k+1);
 end

 eo=mat2str(H);%将矩阵转成字符串
 [h,t]=size(eo);
 eo=eo(h+1:t-1);%去掉开头结尾的[]

% 寻找最大值
function Localmax = Findlocalmax(direction,ix,iy,mag);

[m,n,o] = size(mag);

switch direction
 case 1
 idx = find((iy <= 0 & ix > -iy) | ( iy >= 0 & ix < -iy));
 case 2
 idx = find((ix > 0 & -iy >= ix) | ( ix < 0 & -iy <= ix));
 case 3
 idx = find((ix <= 0 & ix > iy) | ( ix >= 0 & ix < iy));
 case 4
 idx = find((iy < 0 & ix <= iy) | ( iy > 0 & ix >= iy));
end

if ~isempty(idx)
 v = mod(idx,m);
 extIdx = find(v==1|v==0|idx<=m|(idx>(n-1)*m));
 idx(extIdx)=[];
end

ixv = ix(idx);
iyv = iy(idx);
gradmag = mag(idx);
switch direction
 case 1
 d = abs(iyv./ixv);
 gradmag1 = mag(idx + m).*(1-d) + mag(idx + m - 1).*d;
 gradmag2 = mag(idx - m).*(1-d) + mag(idx - m + 1).*d;
 case 2
 d = abs(ixv./iyv);
 gradmag1 = mag(idx - 1).*(1-d) + mag(idx + m - 1).*d;
 gradmag2 = mag(idx + 1).*(1-d) + mag(idx - m + 1).*d;
 case 3
 d = abs(ixv./iyv);
 gradmag1 = mag(idx - 1).*(1-d) + mag(idx - m - 1).*d;
 gradmag2 = mag(idx + 1).*(1-d) + mag(idx + m + 1).*d;
 case 4
 d = abs(iyv./ixv);
 gradmag1 = mag(idx - m).*(1-d) + mag(idx - m - 1).*d;
 gradmag2 = mag(idx + m).*(1-d) + mag(idx + m + 1).*d;
end

Localmax = idx(gradmag >= gradmag1 & gradmag >= gradmag2 );