[MATLAB图像处理] Canny算子提取边缘问题

[xixijida]

[B]谁能帮我调试一下Canny算子提取图像边缘的程序，谢谢

function e=canny_edge(I,sigma)
%function e=edge(I,'canny',thresh,sigma);
%该函数实现Canny算子提取边缘点
%输入图像为I，标准差sigma，输出为边缘图像e
[m,n]=size(I);
Rr=2:m-1;cc=2:n-1;
e=repmat(logical(uint8(0)),m,n);
%产生同样大小的边缘图像e，初始化为1 ，即初始化边缘
GaussianDieOff=-0.001;%设定高斯函数消失门限
PercentOfPixelsNotEdges=-7;%用于计算边缘门限
ThresholdRatio=-4;%设置两个门限的比例
%首先设计高斯滤波器和它的微分
pw=1:30;
%设定滤波器宽度
ssq=sigma*sigma;
%计算方差
width=max(find(exp(-(pw.*pw)/(2*sigma*sigma))>GaussianDieOff));
%计算滤波算子宽度
t=(-width:width);
len=2*width+1;
t3=[t-.5;t;t+.5];
%对每个像素左右各半个像素位置的值进行平均
gau=sum(exp(-(t3.*t3)/(2*ssq))).'/(6*pi*ssq);
%一维高斯滤波器
dgau=(-t.*exp(-(t.*t)/(2*ssq))/ssq).';
%高斯滤波器的微分
ra=size(I,1);
ca=size(I,2);
ay=255*double(I);ax=255*double(I');
h=conv(gau,dgau);
%利用高斯函数滤除噪声和用高斯算子的一阶微分对图像滤波合并为一个算子
ax=conv2(ax,h,'same').';
%产生x方向滤波
ay=conv2(ay,h,'same');
%产生y方向滤波
mag=sqrt((ax.*ax)+(ay.*ay));
%计算滤波结果的幅度
magmax=max(mag();
if magmax>0
mag=mag/magmax;
%对滤波幅度进行归一化
end
%下面根据滤波幅度的概率密度计算滤波门限
[counts,x]=imhist(mag,64);
%计算滤波结果的幅度的直方图
highThresh=min(find(cumsum(counts)>PercentOfPixelsNotEdges*m*n))/64;
%通过设定非边缘点的比例来确定高门限
lowThresh=ThresholdRatio*highThresh;
%设置低门限为高门限乘以比例因子
thresh=[lowThresh,highThresh];
%下面进行非极大抑制
%大于高门限的点归于强边缘图像
%小于低门限的点归于弱边缘图像
idxStrong=[];
for dir=1:4
idxLocalMax=cannyFindLocalMaxima(dir,ax,ay,mag);
idxWeak=idxLocalMax(mag(idxLocalMax)>lowThresh);
e(idxWeak)=1;
idxStrong=[idxStrong;idxWeak(mag(idxWeak)>highThresh)];
end
rstrong=rem(idxStrong-1,m)+1;%rem是求余数
cstrong=floor((idxStrong-1)/m)+1;%向-∞取整
e=bwselect(e,cstrong,rstrong,8);
%通过形态学算子将两幅图像的边缘进行连接
e=bwmorph(e,'thin',1)
%通过形态学算子边缘细化


上面代码中的cannyFindLocalMaxima实现非极大抑制功能，其代码如下：

Function idxLocalMax= cannyFindLocalMaxima(direction,ix,iy,mag);
%输入direction为四个方向,ix为图象在x方向滤波结果.
% iy为图象在y方向滤波结果,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
idxlocalMax=idx(gradmag>=gradmag1&gradmag>=gradmag2);
%进行非极大抑制