圆拟合matlab程序，以及圆心 半径的确定

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;  % 方