通过三个不共线的平面点绘制圆形/MATLAB实现
1: 三个不共线的点求解圆心和半径
假设三个不共线的点的坐标分别为
P
1
(
x
1
,
y
1
)
、
P
2
(
x
2
,
y
2
)
P1(x1,y1)、P2(x2,y2)
P1(x1,y1)、P2(x2,y2)和
P
3
(
x
3
,
y
3
)
P3(x3,y3)
P3(x3,y3),设过三个不共线点的圆的圆心和半径分别为
(
x
,
y
)
和
R
(x,y)和R
(x,y)和R,根据圆的公式可以得到:
{
(
x
−
x
1
)
2
+
(
y
−
y
1
)
2
=
R
2
(
x
−
x
2
)
2
+
(
y
−
y
2
)
2
=
R
2
(
x
−
x
3
)
2
+
(
y
−
y
3
)
2
=
R
2
\begin{cases}(x-x1)^2+(y-y1)^2=R^2\\(x-x2)^2+(y-y2)^2=R^2\\(x-x3)^2+(y-y3)^2=R^2\end{cases}
⎩⎪⎨⎪⎧(x−x1)2+(y−y1)2=R2(x−x2)2+(y−y2)2=R2(x−x3)2+(y−y3)2=R2将上述公式展开移项后得到:
{
x
2
+
y
2
+
x
1
2
+
y
1
2
=
2
x
1
x
+
2
y
1
y
+
R
2
x
2
+
y
2
+
x
2
2
+
y
2
2
=
2
x
2
x
+
2
y
2
y
+
R
2
x
2
+
y
2
+
x
3
2
+
y
3
2
=
2
x
3
x
+
2
y
3
y
+
R
2
\begin{cases}x^2+y^2+x1^2+y1^2=2x1x+2y1y+R^2\\x^2+y^2+x2^2+y2^2=2x2x+2y2y+R^2\\x^2+y^2+x3^2+y3^2=2x3x+2y3y+R^2\end{cases}
⎩⎪⎨⎪⎧x2+y2+x12+y12=2x1x+2y1y+R2x2+y2+x22+y22=2x2x+2y2y+R2x2+y2+x32+y32=2x3x+2y3y+R2通过对上述公式合并整理可以得到:
{
(
x
2
−
x
1
)
x
+
(
y
2
−
y
1
)
y
=
0.5
(
x
2
2
+
y
2
2
−
x
1
2
−
y
1
2
)
(
x
3
−
x
1
)
x
+
(
y
3
−
y
1
)
y
=
0.5
(
x
3
2
+
y
3
2
−
x
1
2
−
y
1
2
)
(
x
3
−
x
2
)
x
+
(
y
3
−
y
2
)
y
=
0.5
(
x
3
2
+
y
3
2
−
x
2
2
−
y
2
2
)
\begin{cases}(x2-x1)x+(y2-y1)y=0.5(x2^2+y2^2-x1^2-y1^2)\\(x3-x1)x+(y3-y1)y=0.5(x3^2+y3^2-x1^2-y1^2)\\(x3-x2)x+(y3-y2)y=0.5(x3^2+y3^2-x2^2-y2^2)\\\end{cases}
⎩⎪⎨⎪⎧(x2−x1)x+(y2−y1)y=0.5(x22+y22−x12−y12)(x3−x1)x+(y3−y1)y=0.5(x32+y32−x12−y12)(x3−x2)x+(y3−y2)y=0.5(x32+y32−x22−y22)将上述方程构成矩阵形式,既有:
A
θ
=
B
A \theta = B
Aθ=B其中
A
=
[
x
2
−
x
1
y
2
−
y
1
x
3
−
x
1
y
3
−
y
1
x
3
−
x
2
y
3
−
y
2
]
,
B
=
0.5
[
x
2
2
+
y
2
2
−
x
1
2
−
y
1
2
x
3
2
+
y
3
2
−
x
1
2
−
y
1
2
x
3
2
+
y
3
2
−
x
2
2
−
y
2
2
]
A=\begin{bmatrix} x2-x1 & y2-y1 \\ x3-x1 & y3-y1 \\ x3-x2 & y3-y2 \end{bmatrix},B=0.5\begin{bmatrix} x2^2+y2^2-x1^2-y1^2 \\ x3^2+y3^2-x1^2-y1^2 \\ x3^2+y3^2-x2^2-y2^2 \end{bmatrix}
A=⎣⎡x2−x1x3−x1x3−x2y2−y1y3−y1y3−y2⎦⎤,B=0.5⎣⎡x22+y22−x12−y12x32+y32−x12−y12x32+y32−x22−y22⎦⎤
由于上述方程为线性方程,也不存在误差,利用最小二乘求得圆心的坐标为:
θ
=
(
A
T
A
)
−
1
A
T
B
\theta =(A^{T}A) ^{-1}A^{T}B
θ=(ATA)−1ATB
那么圆心的坐标有:
x
=
θ
(
1
)
,
y
=
θ
(
2
)
x=\theta(1),y=\theta(2)
x=θ(1),y=θ(2)
利用得到的圆心可以直接求得圆的半径为:
R
=
(
x
−
x
1
)
2
+
(
y
−
y
1
)
2
R=\sqrt{(x-x1)^2+(y-y1)^2}
R=(x−x1)2+(y−y1)2
2: 得到圆的边的坐标
当得到过三点圆的圆心和半径后,可以利用圆的参数方程得到圆的边的坐标,既有: { x = x 0 + R c o s α y = y 0 + R s i n α , 其 中 α ∈ [ 0 , 2 π ] \begin{cases} x=x0+Rcos{\alpha}\\ y=y0+Rsin{\alpha}\end{cases} ,其中 \alpha \in{[0,2\pi]} {x=x0+Rcosαy=y0+Rsinα,其中α∈[0,2π]其中, x x x为圆边上某一点的横坐标, y y y为圆边上某一点的纵坐标, α \alpha α为圆边上某一点对应的角度, x 0 x0 x0和 y 0 y0 y0为该圆的圆心, R R R为该圆的半径。
3: 代码实现-MATLAB
function Result = ThreePoint2Circle(P1, P2, P3)
%% 求圆心和半径
x1 = P1(1); x2 = P2(1); x3 = P3(1);
y1 = P1(2); y2 = P2(2); y3 = P3(2);
z1 = x2^2 + y2^2 - x1^2 - y1^2;
z2 = x3^2 + y3^2 - x1^2 - y1^2;
z3 = x3^2 + y3^2 - x2^2 - y2^2;
A = [(x2-x1), (y2-y1); (x3-x1), (y3-y1); (x3-x2), (y3-y2)];
B = 0.5*[z1; z2; z3];
P0 = (A'*A)\A'*B;
R1 = sqrt( (P0(1) - P1(1))^2 + (P0(2) - P1(2))^2 );
R2 = sqrt( (P0(1) - P2(1))^2 + (P0(2) - P2(2))^2 );
R3 = sqrt( (P0(1) - P3(1))^2 + (P0(2) - P3(2))^2 );
R = (R1 + R2 + R3)/3;
%% 绘制圆
theta = (0:pi/360:2*pi)';
Result = zeros(size(theta,1),4);
for i = 1: size(theta,1)
Result(i,1) = i;
Result(i,2) = theta(i);
Result(i,3) = P0(1) + R*cos(theta(i));
Result(i,4) = P0(2) + R*sin(theta(i));
end
figure();plot(Result(:,3),Result(:,4),'b-');hold on;
grid on; xlabel('x');ylabel('y'); axis equal;
end
输出的结果为0.5°间隔,其中第1列为点号,第2列为角度,第3列为横坐标,第4列为纵坐标。
4: 代码测试
(1)过 ( 1 , 1 ) 、 ( 3 , 2 ) 、 ( 5 , 1 ) (1,1)、(3,2)、(5,1) (1,1)、(3,2)、(5,1)画圆
(2)过 ( 1.2 , 1.2 ) 、 ( 8 , 0.5 ) 、 ( 13.8 , 1.2 ) (1.2,1.2)、(8,0.5)、(13.8,1.2) (1.2,1.2)、(8,0.5)、(13.8,1.2)画圆