拟合圆并求圆心(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;