最小二乘圆拟合模型公式推导&MATLAB代码求圆拟合

理论部分1：
https://wenku.baidu.com/view/ecda32525beef8c75fbfc77da26925c52dc59156.html
最小二乘圆拟合模型公式推导
理论部分2：
https://blog.csdn.net/liyuanbhu/article/details/52891868
三点确定一个圆的计算方法
理论部分3：
https://blog.csdn.net/Sppy_z/article/details/104877864?utm_medium=distribute.pc_relevant.none-task-blog-2_defaultbaidujs_title~default-0-104877864-blog-6607930.pc_relevant_antiscanv3&spm=1001.2101.3001.4242.1&utm_relevant_index=3
由圆上三点确定圆心和半径

代码部分：采用一般方程直接拟合
https://wenku.baidu.com/view/f38f72f36c1aff00bed5b9f3f90f76c661374c94.html
MATLAB代码求圆拟合

代码部分2：采用公式推导拟合
https://www.docin.com/p-260086637.html
matlab代码公式推导拟合

稍作整理，直接拟合：

%%
% https://wenku.baidu.com/view/ecda32525beef8c75fbfc77da26925c52dc59156.html
% https://wenku.baidu.com/view/f38f72f36c1aff00bed5b9f3f90f76c661374c94.html
%%
clear;close all;clc;
format compact;
figure
t=0:0.01:pi;
a=20;
b=30;
r=5;
x=a+r*cos(t)+randn(1,315);
y=b+r*sin(t)+randn(1,315);
plot(x,y,'b.')
hold on;grid on;
x=x(:); % Transpose
y=y(:); 
Mole=[x y ones(size(x))];
Deno=[-(x.^2+y.^2)];
Result1=Mole\Deno
% x^2+y^2+ax+by+c=0 ---> ax+by+c=-(x^2+y^2)
% [x y 1]*[a;b;c]=-(x^2+y^2)
% [a;b;c]=[x y 1]\[-(x^2+y^2)]
% Column principal elimination method
% Least square fitting can also be used
m=[x y ones(size(x))]\[-(x.^2+y.^2)] 
xc=-0.5*m(1);
yc=-0.5*m(2);
R=sqrt((m(1)^2+m(2)^2)/4-m(3))
plot(xc,yc,'r-x',(xc+R*cos(t)),(yc+R*sin(t)),'r.');
axis equal;

稍作整理，公式推导拟合：

%% https://www.docin.com/p-260086637.html
%%
function [R,xc,yc]=circle(x,y,N)

clc;close all;clear;
t=0:0.01:pi;
a=20;b=30;r=5;
% point(x,y)
x=a+r*cos(t)+randn(1,315);
y=b+r*sin(t)+randn(1,315);
figure(1)
plot(x,y,'b.')
hold on;grid on;axis equal;
% (x-xc)^2+(y-yc)^2=R^2
% x^2+y^2+ax+by+c=0
% xc=-a/2  yc=-b/2  R=sqrt(a^2+b^2-4c)/2
N=315; % Number of points
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
 % Center and radius
 xc=a/(-2)
 yc=b/(-2)
 R=sqrt(a*a+b*b-4*c)/2

 plot(xc,yc,'r-x',(xc+R*cos(t)),(yc+R*sin(t)),'r.');


end