1.圓方程定義
常規(guī)圓方程定義為
(
x
?
x
0
)
2
+
(
y
?
y
0
)
2
=
r
2
(x-x_0)^2+(y-y_0)^2 = r^2
(x?x0?)2+(y?y0?)2=r2
可以改寫為
x
2
+
y
2
+
a
x
+
b
y
+
c
=
0.
x^2+y^2+ax+by+c=0.
x2+y2+ax+by+c=0.
其中
a
=
?
2
x
0
,
b
=
?
2
y
0
,
c
=
x
0
2
+
y
0
2
?
r
2
a=-2x_0,b=-2y_0,c={x_0}^2+{y_0}^2-r^2
a=?2x0?,b=?2y0?,c=x0?2+y0?2?r2文章來(lái)源：http://www.zghlxwxcb.cn/news/detail-421271.html
2.最小能量函數(shù)定義
要使得圓方程最為準(zhǔn)確，則要是所有的點(diǎn)盡可能滿足方程，因此需使函數(shù)F
F
=
∑
(
x
i
2
+
y
i
2
+
a
x
i
+
b
y
i
+
c
)
2
=
0
,
i
∈
[
1
,
N
]
,
N
為
點(diǎn)
數(shù)
目
F=\sum( {x_i}^2+{y_i}^2+a{x_i}+b{y_i}+c)^2 =0 , i\in [1,N],N為點(diǎn)數(shù)目
F=∑(xi?2+yi?2+axi?+byi?+c)2=0,i∈[1,N],N為點(diǎn)數(shù)目
但實(shí)際上其不可能為0，因此根據(jù)最小二乘法求解其最小值
對(duì)函數(shù)F求導(dǎo)可得
?
F
a
=
∑
2
(
x
i
2
+
y
i
2
+
a
x
i
+
b
y
i
+
c
)
?
x
i
\frac{\partial{F}}{a} = \sum{2( {x_i}^2+{y_i}^2+a{x_i}+b{y_i}+c)*x_i}
a?F?=∑2(xi?2+yi?2+axi?+byi?+c)?xi?
?
F
b
=
∑
2
(
x
i
2
+
y
i
2
+
a
x
i
+
b
y
i
+
c
)
?
y
i
\frac{\partial{F}} = \sum{2( {x_i}^2+{y_i}^2+a{x_i}+b{y_i}+c)*y_i}
b?F?=∑2(xi?2+yi?2+axi?+byi?+c)?yi?
?
F
c
=
∑
2
(
x
i
2
+
y
i
2
+
a
x
i
+
b
y
i
+
c
)
\frac{\partial{F}}{c} = \sum{2( {x_i}^2+{y_i}^2+a{x_i}+b{y_i}+c)}
c?F?=∑2(xi?2+yi?2+axi?+byi?+c)
當(dāng)這些導(dǎo)數(shù)為0是，函數(shù)最小
為方便求解，將其寫成a,b,c的矩陣函數(shù)形式，則可得
接下來(lái)便是求解方程AX=B文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-421271.html
3.matlab代碼
- 實(shí)例
% 圓擬合實(shí)例
clc;
close all;
%% 創(chuàng)建測(cè)試數(shù)據(jù)
num =1000;
% 噪聲
rng(1);
noise1=0.2*randn(num,1);
rng(2);
noise2=0.2*randn(num,1);
% 創(chuàng)建初始數(shù)據(jù)
theta=(2*rand(num,1)-1)*pi;
center = [10,10];
r=3;
x= center(1) + r*cos(theta)+noise1;
y= center(2) + r*sin(theta)+noise2;
%% 擬合圓
figure;
scatter(x,y);
hold on;
[a,b,c] = fit_circle(x,y);
center_cal=[-a/2,-b/2];
r_cal = sqrt(a*a/4+b*b/4-c);
theta_cal =0:0.01:2*pi;
x_cal= center_cal(1) + r_cal*cos(theta_cal);
y_cal= center_cal(2) + r_cal*sin(theta_cal);
plot(x_cal,y_cal);
- 圓擬合函數(shù)
function [a,b,c] = fit_circle(x,y)
% function fit_circle.m
% -------------------------------------------------------------------------
% 圓擬合
% 圓方程按照x^2+y^2+ax+by+c=0
% 最小化函數(shù)定義為 sum((x_i)^2+(y_i)^2+a(x_i)+b(y_i)+c),i為1到n的整數(shù)
% 計(jì)算導(dǎo)數(shù)，化為AX=B的形式，X為[a;b;c]
% Inputs:
% x = 點(diǎn)x坐標(biāo),寫成(1000,1)豎排格式
% y = 點(diǎn)y坐標(biāo)
% Output:
% a,b,c = 圓方程對(duì)應(yīng)參數(shù)
% Last updated: November 14, 2022
% 構(gòu)建A
A1_1 = sum(x.^2);
A1_2 = sum(x.*y);
A1_3 = sum(x);
A2_1 = A1_2;
A2_2 = sum(y.^2);
A2_3 = sum(y);
A3_1 = A1_3;
A3_2 = A2_3;
A3_3 = size(x,1);
A=[A1_1,A1_2,A1_3;
A2_1,A2_2,A2_3;
A3_1,A3_2,A3_3];
% 構(gòu)建B
B1 = -sum(x.^3+(y.^2).*x);
B2 = -sum(y.^3+(x.^2).*y);
B3 = -sum(x.^2+y.^2);
B = [B1;B2;B3];
% 計(jì)算
X = A\B;
a=X(1);
b=X(2);
c=X(3);
end
到了這里，關(guān)于標(biāo)準(zhǔn)二維圓擬合(matlab代碼)的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
