一、问题描述
给定三维空间不共面的四个点 A , B , C , D A,B,C,D A,B,C,D，求由此四点确定的球面方程 ( x − x 0 ) 2 + ( y − y 0 ) 2 + ( z − z 0 ) 2 = R 2 (x-x_0)^2+(y-y_0)^2+(z-z_0)^2=R^2 (x−x0)2+(y−y0)2+(z−z0)2=R2。
二、推导步骤
方法一
设坐标
A
(
x
1
,
y
1
,
z
1
)
,
B
(
x
2
,
y
2
,
z
2
)
,
C
(
x
3
,
y
3
,
z
3
)
,
D
(
x
4
,
y
4
,
z
4
)
A(x_1,y_1,z_1),B(x_2,y_2,z_2),C(x_3,y_3,z_3),D(x_4,y_4,z_4)
A(x1,y1,z1),B(x2,y2,z2),C(x3,y3,z3),D(x4,y4,z4)，由于此四点在球面上，因而：
{
(
x
1
−
x
0
)
2
+
(
y
1
−
y
0
)
2
+
(
z
1
−
z
0
)
2
=
R
2
(
1.1
)
(
x
2
−
x
0
)
2
+
(
y
2
−
y
0
)
2
+
(
z
2
−
z
0
)
2
=
R
2
(
1.2
)
(
x
3
−
x
0
)
2
+
(
y
3
−
y
0
)
2
+
(
z
3
−
z
0
)
2
=
R
2
(
1.3
)
(
x
4
−
x
0
)
2
+
(
y
4
−
y
0
)
2
+
(
z
4
−
z
0
)
2
=
R
2
(
1.4
)
(1)
\begin{cases} (x_1-x_0)^2+(y_1-y_0)^2+(z_1-z_0)^2=R^2 & (1.1)\\ (x_2-x_0)^2+(y_2-y_0)^2+(z_2-z_0)^2=R^2 & (1.2)\\ (x_3-x_0)^2+(y_3-y_0)^2+(z_3-z_0)^2=R^2 & (1.3)\\ (x_4-x_0)^2+(y_4-y_0)^2+(z_4-z_0)^2=R^2 & (1.4)\\ \tag 1 \end{cases}
⎩
⎨
⎧(x1−x0)2+(y1−y0)2+(z1−z0)2=R2(x2−x0)2+(y2−y0)2+(z2−z0)2=R2(x3−x0)2+(y3−y0)2+(z3−z0)2=R2(x4−x0)2+(y4−y0)2+(z4−z0)2=R2(1.1)(1.2)(1.3)(1.4)(1)
式
(
1.1
)
(1.1)
(1.1)至
(
1.4
)
(1.4)
(1.4)两两做差，整理可得：
(
x
j
−
x
i
)
x
0
+
(
y
j
−
y
i
)
y
0
+
(
z
j
−
z
i
)
z
0
=
[
(
x
j
−
x
i
)
(
x
j
+
x
i
)
+
(
y
j
−
y
i
)
(
y
j
+
y
i
)
+
(
z
j
−
z
i
)
(
z
j
+
z
i
)
]
/
2
(2)
\begin{aligned} & (x_j - x_i)x_0 + (y_j - y_i)y_0 + (z_j - z_i)z_0 \\ &= [(x_j-x_i)(x_j+x_i)+(y_j-y_i)(y_j+y_i) \\ & +(z_j-z_i)(z_j+z_i)]/2 \tag 2 \\ \end{aligned}
(xj−xi)x0+(yj−yi)y0+(zj−zi)z0=[(xj−xi)(xj+xi)+(yj−yi)(yj+yi)+(zj−zi)(zj+zi)]/2(2)
其中，
i
,
j
=
1
,
2
,
3
,
4
i,j=1,2,3,4
i,j=1,2,3,4且
i
≠
j
i\neq j
i=j，即式
(
2
)
(2)
(2)有
C
4
2
=
6
C_4^2=6
C42=6个方程。
6个方程中任取3个组成方程组，共有
C
6
3
=
20
C_6^3=20
C63=20种情况，不妨选取以下方程组求解球心：
M
p
=
N
(3)
Mp=N \tag 3
Mp=N(3)
其中，
M
=
[
x
2
−
x
1
y
2
−
y
1
z
2
−
z
1
x
3
−
x
1
y
3
−
y
1
z
4
−
z
1
x
4
−
x
1
y
4
−
y
1
z
4
−
z
1
]
(4)
M=\left[ \begin{matrix} x_2 - x_1 & y_2 - y_1 & z_2 - z_1 \\ x_3 - x_1 & y_3 - y_1 & z_4 - z_1 \\ x_4 - x_1 & y_4 - y_1 & z_4 - z_1 \\ \end{matrix} \right] \tag 4
M=
x2−x1x3−x1x4−x1y2−y1y3−y1y4−y1z2−z1z4−z1z4−z1
(4)
N
=
[
[
(
x
2
−
x
1
)
(
x
2
+
x
1
)
+
(
y
2
−
y
1
)
(
y
2
+
y
1
)
+
(
z
2
−
z
1
)
(
z
2
+
z
1
)
]
/
2
[
(
x
3
−
x
1
)
(
x
3
+
x
1
)
+
(
y
3
−
y
1
)
(
y
3
+
y
1
)
+
(
z
3
−
z
1
)
(
z
3
+
z
1
)
]
/
2
[
(
x
4
−
x
1
)
(
x
4
+
x
1
)
+
(
y
4
−
y
1
)
(
y
4
+
y
1
)
+
(
z
4
−
z
1
)
(
z
4
+
z
1
)
]
/
2
]
(5)
N=\left[ \begin{matrix} [(x_2-x_1)(x_2+x_1)+(y_2-y_1)(y_2+y_1) +(z_2-z_1)(z_2+z_1)]/2 \\ [(x_3-x_1)(x_3+x_1)+(y_3-y_1)(y_3+y_1) +(z_3-z_1)(z_3+z_1)]/2 \\ [(x_4-x_1)(x_4+x_1)+(y_4-y_1)(y_4+y_1) +(z_4-z_1)(z_4+z_1)]/2 \\ \end{matrix} \right] \tag 5
N=
[(x2−x1)(x2+x1)+(y2−y1)(y2+y1)+(z2−z1)(z2+z1)]/2[(x3−x1)(x3+x1)+(y3−y1)(y3+y1)+(z3−z1)(z3+z1)]/2[(x4−x1)(x4+x1)+(y4−y1)(y4+y1)+(z4−z1)(z4+z1)]/2
(5)
经计算，矩阵
M
M
M的行列式值满足：
∣
M
∣
=
(
[
x
2
−
x
1
y
2
−
y
1
z
2
−
z
1
]
×
[
x
3
−
x
1
y
3
−
y
1
z
3
−
z
1
]
)
⋅
[
x
4
−
x
1
y
4
−
y
1
z
4
−
z
1
]
=
(
A
B
⃗
×
A
C
⃗
)
⋅
A
D
⃗
(6)
\begin{aligned} |M|&=(\left[ \begin{matrix} x_2 - x_1\\ y_2 - y_1 \\ z_2 - z_1 \\ \end{matrix} \right]\times \left[ \begin{matrix} x_3 - x_1\\ y_3 - y_1 \\ z_3 - z_1 \\ \end{matrix} \right])\cdot \left[ \begin{matrix} x_4 - x_1\\ y_4 - y_1 \\ z_4 - z_1 \\ \end{matrix} \right] \\ & =(\vec{AB}\times\vec{AC})\cdot\vec{AD} \end{aligned} \tag 6
∣M∣=(
x2−x1y2−y1z2−z1
×
x3−x1y3−y1z3−z1
)⋅
x4−x1y4−y1z4−z1
=(AB×AC)⋅AD(6)
当
A
,
B
,
C
,
D
A,B,C,D
A,B,C,D不共面时，混合积
(
A
B
⃗
×
A
C
⃗
)
⋅
A
D
⃗
(\vec{AB}\times\vec{AC})\cdot\vec{AD}
(AB×AC)⋅AD为四面体
A
B
C
D
ABCD
ABCD的有向体积，不为零，因此矩阵
M
M
M可逆，矩阵方程
(
3
)
(3)
(3)有唯一解。
方法二
将球面方程展开成：
x
2
+
y
2
+
z
2
−
2
x
0
x
−
2
y
0
y
−
2
z
0
z
+
x
0
2
+
y
0
2
+
z
0
2
−
R
2
=
0
(7)
x^2+y^2+z^2-2x_0x-2y_0y-2z_0z+x_0^2+ y_0^2+ z_0^2-R^2=0 \tag 7
x2+y2+z2−2x0x−2y0y−2z0z+x02+y02+z02−R2=0(7)
根据多项式的项和系数特征可以将式(7)写成：
x
2
+
y
2
+
z
2
+
a
x
+
b
y
+
c
z
+
d
=
0
(8)
x^2+y^2+z^2+ax+by+cz+d=0 \tag 8
x2+y2+z2+ax+by+cz+d=0(8)
将给定4个点
(
x
i
,
y
i
,
z
i
)
(
i
=
1
,
2
,
3
,
4
)
(x_i,y_i,z_i)(i=1,2,3,4)
(xi,yi,zi)(i=1,2,3,4)带入式(8)，可以得到关于待定系数
a
,
b
,
c
,
d
a,b,c,d
a,b,c,d的四元一次线性方程组。求解四元一次线性方程组，可以求得球心坐标：
{
x
0
=
−
a
/
2
y
0
=
−
b
/
2
z
0
=
−
c
/
2
(9)
\begin{cases} x_0=-a/2 &\\ y_0=-b/2 &\\ z_0=-c/2 &\\ \tag 9 \end{cases}
⎩
⎨
⎧x0=−a/2y0=−b/2z0=−c/2(9)
球半径：
R
=
(
a
/
2
)
2
+
(
b
/
2
)
2
+
(
c
/
2
)
2
−
d
(10)
R=\sqrt{(a/2)^2+(b/2)^2+(c/2)^2-d} \tag {10}
R=(a/2)2+(b/2)2+(c/2)2−d(10)
三、MATLAB代码
%{
Function: solve_sphere_params
Description: 求四点确定的球面参数
Input: 空间不共面四个点A,B,C,D
Output: 球面球心sphereCenter, 半径radius, 求解状态status(1表示成功，0表示失败)
Author: Marc Pony(marc_pony@163.com)
%}
function [sphereCenter, radius, status] = solve_sphere_params(A, B, C, D)
x1 = A(1);
y1 = A(2);
z1 = A(3);
x2 = B(1);
y2 = B(2);
z2 = B(3);
x3 = C(1);
y3 = C(2);
z3 = C(3);
x4 = D(1);
y4 = D(2);
z4 = D(3);
sphereCenter = zeros(3, 1);
radius = 0.0;
status = 1;
if abs((y1 - y4) * ((x1 - x2) * (z1 - z3) - (x1 - x3) * (z1 - z2)) ...
- (z1 - z4) * ((x1 - x2) * (y1 - y3) - (x1 - x3) * (y1 - y2)) ...
- (x1 - x4) * ((y1 - y2) * (z1 - z3) - (y1 - y3) * (z1 - z2))) < 1.0e-8 %dot(cross(AB, AC), AD)
status = 0;
return;
else
a11 = x2 - x1;
a12 = y2 - y1;
a13 = z2 - z1;
b1 = 0.5 * ((x2 - x1) * (x2 + x1) + (y2 - y1) * (y2 + y1) + (z2 - z1) * (z2 + z1));
a21 = x3 - x1;
a22 = y3 - y1;
a23 = z3 - z1;
b2 = 0.5 * ((x3 - x1) * (x3 + x1) + (y3 - y1) * (y3 + y1) + (z3 - z1) * (z3 + z1));
a31 = x4 - x1;
a32 = y4 - y1;
a33 = z4 - z1;
b3 = 0.5 * ((x4 - x1) * (x4 + x1) + (y4 - y1) * (y4 + y1) + (z4 - z1) * (z4 + z1));
temp = a11 * (a22 * a33 - a23 * a32) + a12 * (a23 * a31 - a21 * a33) + a13 * (a21 * a32 - a22 * a31);
x0 = ((a12 * a23 - a13 * a22) * b3 + (a13 * a32 - a12 * a33) * b2 + (a22 * a33 - a23 * a32) * b1) / temp;
y0 = -((a11 * a23 - a13 * a21) * b3 + (a13 * a31 - a11 * a33) * b2 + (a21 * a33 - a23 * a31) * b1) / temp;
z0 = ((a11 * a22 - a12 * a21) * b3 + (a12 * a31 - a11 * a32) * b2 + (a21 * a32 - a22 * a31) * b1) / temp;
radius = sqrt((x0 - x1)^2 + (y0 - y1)^2 + (z0 - z1)^2);
sphereCenter = [x0; y0; z0];
end
end
%{
Function: solve_sphere_params2
Description: 求四点确定的球面参数
Input: 空间不共面四个点A,B,C,D
Output: 球面球心sphereCenter, 半径radius, 求解状态status(1表示成功，0表示失败)
Author: Marc Pony(marc_pony@163.com)
%}
function [sphereCenter, radius, status] = solve_sphere_params2(A, B, C, D)
x1 = A(1);
y1 = A(2);
z1 = A(3);
x2 = B(1);
y2 = B(2);
z2 = B(3);
x3 = C(1);
y3 = C(2);
z3 = C(3);
x4 = D(1);
y4 = D(2);
z4 = D(3);
A = [x1, y1, z1, 1; ...
x2, y2, z2, 1; ...
x3, y3, z3, 1; ...
x4, y4, z4, 1
];
B = -[x1^2 + y1^2 + z1^2; x2^2 + y2^2 + z2^2; x3^2 + y3^2 + z3^2; x4^2 + y4^2 + z4^2];
if abs(det(A)) > 1.0e-6
X = A \ B;
a = X(1);
b = X(2);
c = X(3);
d = X(4);
sphereCenter = [-a/2; -b/2; -c/2];
radius = sqrt((a/2)^2 + (b/2)^2 + (c/2)^2 - d);
status = 1;
else
status = 0;
sphereCenter = [];
radius = [];
return;
end
end
%{
Function: draw_sphere
Description: 画球面
Input: 球心sphereCenter，球半径radius
Output: 无
Author: Marc Pony(marc_pony@163.com)
%}
function draw_sphere(sphereCenter, radius)
[x,y,z] = sphere(200);
x = x * radius + sphereCenter(1);
y = y * radius + sphereCenter(2);
z = z * radius + sphereCenter(3);
h = surf(x, y, z);
xlabel('x')
ylabel('y')
zlabel('z')
set(h, 'FaceAlpha', 0.3, 'MarkerEdgeColor', 'none')
shading interp
end
clc
clear
close all
pos = 10.0 * rand(3, 4);
A = pos(:, 1);
B = pos(:, 2);
C = pos(:, 3);
D = pos(:, 4);
[sphereCenter, radius, status] = solve_sphere_params(A, B, C, D)
[sphereCenter2, radius2, status2] = solve_sphere_params2(A, B, C, D)
err = [abs( (A(1) - sphereCenter(1))^2 + (A(2) - sphereCenter(2))^2 + (A(3) - sphereCenter(3))^2 - radius^2 ); ...
abs( (B(1) - sphereCenter(1))^2 + (B(2) - sphereCenter(2))^2 + (B(3) - sphereCenter(3))^2 - radius^2 ); ...
abs( (C(1) - sphereCenter(1))^2 + (C(2) - sphereCenter(2))^2 + (C(3) - sphereCenter(3))^2 - radius^2 ); ...
abs( (D(1) - sphereCenter(1))^2 + (D(2) - sphereCenter(2))^2 + (D(3) - sphereCenter(3))^2 - radius^2 )]
figure('color', 'w')
draw_sphere(sphereCenter, radius)
hold on
plot3(pos(1, :), pos(2, :), pos(3, :), 'r+')
axis equal tight
扫一扫
1679
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
