最小二乘法擬合圓

方法並非原創，網上找的資料參考

通過編寫代碼熟悉Matlab算法仿真。

今后需要改進的地方：

1. 克服噪聲的方法
2. 實時數據的選擇
3. 解決實際應用問題

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);
    B = b/(-2);
    R = sqrt(a*a + b*b - 4*c)/2;

end

% Main Test Code

clear all;
np = 20;
noise = 2;
angle = (rand(1,np)-0.5)*2*pi;
a = -22;     %center of circle X 
b = 33;     %center of circle Y
r = 25;      %radius
x = a+r*cos(angle)+randn(1,np)*noise;    %random matrix(1x315) becuase pi=3.14
y = b+r*sin(angle)+randn(1,np)*noise;
plot(x,y,'g*');
hold on;
[R,xc,yc] = circ(x,y,np);
plot(xc,yc,'b+',(xc+R*cos(angle)),(yc+R*sin(angle)),'r*');
axis equal;
grid on;