拟合圆并求圆心(matlab)

标签：

clc
 clear all
 I=imread(‘a003.bmp‘);
 [m,n,p]=size(I);
 I=im2bw(I);
 se=strel(‘square‘,3);
 Ia=imerode(I,se);
 Iout=I-Ia;
 figure(1);
 imshow(Iout);
 Iout1=Iout;
 N=1;
 for i=1:m
     for j=1:n
         if Iout(i,j)==1
             p(N,1)=i;
             p(N,2)=j;
             N=N+1;
             Iout(i,j+5:end)=0;
         end
        
     end
 end
 figure(2);
 imshow(Iout);
 for i=1:m
     for j=n:(-1):1
         if Iout1(i,j)==1
             p(N,1)=i;
             p(N,2)=j;
             N=N+1;
             Iout1(i,1:j-5)=0;
         end
     end
 end
 figure(3);
 imshow(Iout1);
 x=p(:,1);
 y=p(:,2);
 [R,x0,y0]=circ(x,y,N-1);
 disp([‘这个圆的半径是: ‘ num2str(R),‘ mm‘ ]);
 disp([‘这个圆圆心坐标:( ‘ num2str(x0) ,‘,‘,num2str(y0),‘)‘]);
 figure(4)
 imshow(I);
 hold on
 plot(x0,y0,‘r+‘);

%%%%%%%%%%%%%%%%%%%%%%%

function [R,A,B]=circ(x,y,N)
x1 = 0;
x2 = 0;
x3 = 0;
y1 = 0;
y2 = 0;
y3 = 0;
x1y1 = 0;
x1y2 = 0;
x2y1 = 0;
for i = 1 : N
x1 = x1 + x(i);
x2 = x2 + x(i)*x(i);
x3 = x3 + x(i)*x(i)*x(i);
y1 = y1 + y(i);
y2 = y2 + y(i)*y(i);
y3 = y3 + y(i)*y(i)*y(i);
x1y1 = x1y1 + x(i)*y(i);
x1y2 = x1y2 + x(i)*y(i)*y(i);
x2y1 = x2y1 + x(i)*x(i)*y(i);
end
C = N * x2 - x1 * x1;
D = N * x1y1 - x1 * y1;
E = N * x3 + N * x1y2 - (x2 + y2) * x1;
G = N * y2 - y1 * y1;
H = N * x2y1 + N * y3 - (x2 + y2) * y1;
a = (H * D - E * G)/(C * G - D * D);
b = (H * C - E * D)/(D * D - G * C);
c = -(a * x1 + b * y1 + x2 + y2)/N;
A = a/(-2); %x 坐标
B = b/(-2); %y 坐标
R = sqrt(a * a + b * b - 4 * c)/2;

拟合圆并求圆心(matlab)

标签：

原文地址：http://blog.csdn.net/u012374012/article/details/44828059