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);
%通過形態學算子將兩幅圖像的邊緣進行連接