%本程序用于實驗canny濾波算法
%作者：大指揮官
%程序狀態：完結，但是效果不如MATLAB自帶的Canny濾波器效果好。

pic=imread（‘e:\test-pic\test3.jpg‘）;
pic=rgb2gray（pic）;
%pic=imresize（pic0.5）;
pic_s=pic;
px=size（pic1）;
py=size（pic2）;
%?figure（1）;
%?imshow（pic）;
%?title（‘原圖‘）;
%%%%%%%%%%%%%%%%%%
%構建高斯濾波算子
%%%%%%%%%%%%%%%%%%
sigma=1;
Gx=floor（3*sigma+0.5）;
x=-1*Gx:1*Gx;
y=x;
[xy]=meshgrid（xy）;
z=1/（2*pi*sigma^2）*exp（-1*（x.^2+y.^2）/（2*sigma.^2））;
????%高斯算子歸一化
z_sum=sum（sum（z））;
z=z/z_sum;
%?figure（2）;
%?mesh（z）;
pic=uint8（conv2（double（pic）z））;
pic=pic（Gx+1:px+GxGx+1:py+Gx）;
pic=uint8（pic）;
figure（3）;
imshow（pic）;
title（strcat（‘sigma=‘num2str（sigma）））;

%%%%%%%%%%%%%%%%%%%%%%%
%計算兩個方向上的一階差分
%%%%%%%%%%%%%%%%%%%%%%%
A=-1/2*[-1?1;-1?1];
B=-1/2*[1?1;-1?-1];
pic_x=conv2（double（pic）A）;
pic_x=pic_x（2:px+12:py+1）;
pic_y=conv2（double（pic）B）;
pic_y=pic_y（2:px+12:py+1）;
%?figure（4）;
%?imshow（pic_x[]）;
%?title（‘x方向差分‘）;
%?figure（5）;
%?imshow（pic_y[]）;
%?title（‘y方向差分‘）;

%%%%%%%%%%%%%%%%%%%%%%%
%計算梯度的幅值與方向
%%%%%%%%%%%%%%%%%%%%%%%
gradamp=sqrt（double（pic_x.^2）+double（pic_y.^2））;
gradori=atan（double（pic_y）./double（pic_x））;
figure（6）;
imshow（gradamp[]）;
title（‘梯度幅值‘）;

%%%%%%%%%%%%%%%%%%%%%%%
%非極大值抑制
%%%%%%%%%%%%%%%%%%%%%%%
%此處是以梯度圖為模板，比較梯度圖在梯度方向上的局部最大值，從而細化邊緣
picre=zeros（size（gradamp））;
picre（1:）=0;
picre（px:）=0;
picre（:1）=0;
picre（:py）=0;
%g1=0;g2=0;g3=0;g4=0;%這四個變量用于插值計算
t1=0;t2=0;??????????%這兩個變量保存插好的值
h1=0;h2=0;??????????%插值計算時的兩個端點的權重
for?i=2:px-1
????for?j=2:py-1
????????if?gradamp（ij）~=0
????????????theta=gradori（ij）;
????????????if?theta<0
????????????????theta=theta+pi;
????????????end
????????????%??g1????g2
????????????%????????c
????????????%????????g3????g4
????????????if?（theta>=pi/2）&&（theta????????????????h1=tan（theta-pi/2）;
????????????????h2=1-h1;
????????????????t1=h1*gradamp（i-1j-1）+h2*gradamp（i-1j）;
????????????????t2=h1*gradamp（i+1j+1）+h2*gradamp（i+1j）;
????????????end
????????????%??g1
????????????%??g2????c????g3
????????????%?????????????g4
????????????if?（theta>=pi*3/4）&&（theta<=pi）
????????????????h1=tan（pi-theta）;
????????????????h2=1-h1;
????????????????t1=h1*g