[replyview]
function e=canny_edge(I,sigma)
%functione=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);     
%通過形態學算子將兩幅圖像的邊緣進行連接