function?canny
????????clear;clc
????????I?=?imread（‘16576.jpg‘）;
????????C=rgb2gray（I）;
????????se=strel（‘disk‘17）;
????????K=imerode（Cse）;
????????O=medfilt2（K[10?10]）;
????????figure（5）
????????imshow（K）
????????BW?=?canny_edge（O10）;
????????figure（7）
????????imshow（BW）
????????r=I（::1）;
????????g=I（::2）;
????????b=I（::3）;
????????
????????bw=imfill（BW‘hole‘）;
????????%figure
????????%imshow（bw）
????????se=strel（‘disk‘4）;
????????bw=bwmorph（bw‘spur‘）;
????????bw=bwmorph（bw‘clean‘inf）;
????????
????????bw=imdilate（bwse）;
????????bw=imfill（bw‘hole‘）;
????????
????????%figure（4）
????????%imshow（bw）
????????bw=imerode（bwse）;
????????L=bwperim（bw4）;
????????r（L）=0;
????????g（L）=255;
????????b（L）=0;
????????figure（3）
????????imshow（cat（3rgb））;
????????hold?on
????????
????????
????????[xy]=find（L==1）;
????????D=[xy];
????????[xy]
????????k=sort（y）;
????????m=size（x1）;
????????n=size（y1）;
????????temp=k（n-901）;
????????j1=0;
????????for?i=1:n
????????????if（y（i1）>=temp）
????????????????D（i:）=0;
????????????????j1=j1+1;
????????????????n=n-1;
????????????else
????????????????continue
????????????end
????????end
????????x=D（:1）;
????????y=D（:2）;
????????x1=sum（x）/n;
????????x2=sum（x‘*x）/n;
????????x3=sum（x‘*diag（x*x‘））/n;
????????
????????y1=sum（y）/n;
????????y2=sum（y‘*y）/n;
????????y3=sum（y‘*diag（y*y‘））/n;
????????
????????xy=sum（x‘*y）/n;
????????xy1=sum（y‘*diag（x*x‘））/n;
????????xy2=sum（x‘*diag（y*y‘））/n;
????????
????????h=y1*y1-y2;
????????i=x1*x1-x2;
????????j=x1*y1-xy;
????????
????????u=x2*x1+x1*y2-x3-xy2;
????????v=y2*y1+y1*x2-y3-xy1;
????????w=i*h+j*j;
????????
???????
????????b=（u*h+v*j）/（2*w）;
????????a=（v*i-u*j）/（2*w）;
????????r=sqrt（b*b-2*x1*b+x2+a*a-2*y1*a+y2）;
????????t=linspace（-pipi）;
????????x=sin（t）‘*r+a;
????????y=cos（t）‘*r+b;
????????a
????????b
????????r
????????
????????plot（xy‘w-‘）;
????????plot（ab‘o-‘）
????????%axis?on
????????
????????
????????%figure（2）
????????%imshow（BW）
????????

%?Canny邊緣檢測的函數
%?Input:
%???a:?input?image
%???sigma:?Gaussian的均方差
function?e=canny_edge（asigma）

a?=?double（a）;?????%?將圖像像素數據轉換為double類型
[mn]?=?size（a）;
e?=?repmat（logical（uint8（0））mn）;??%?生成初始矩陣

OffGate?=?0.0001;??
%Perc?=?0.7;
Th?=?1.1;

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