一、问题描述

给定三维空间不共面的四个点A,B,C,DA,B,C,DA,B,C,D，求由此四点确定的球面方程(x−x0)2+(y−y0)2+(z−z0)2=R2(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(x1,y1,z1),B(x2,y2,z2),C(x3,y3,z3),D(x4,y4,z4)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)，由于此四点在球面上，因而：
{(x1−x0)2+(y1−y0)2+(z1−z0)2=R2(1.1)(x2−x0)2+(y2−y0)2+(z2−z0)2=R2(1.2)(x3−x0)2+(y3−y0)2+(z3−z0)2=R2(1.3)(x4−x0)2+(y4−y0)2+(z4−z0)2=R2(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)两两做差，整理可得：
(xj−xi)x0+(yj−yi)y0+(zj−zi)z0=[(xj−xi)(xj+xi)+(yj−yi)(yj+yi)+(zj−zi)(zj+zi)]/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,4i,j=1,2,3,4i,j=1,2,3,4且i≠ji\neq ji=j，即式(2)(2)(2)有C42=6C_4^2=6C42=6个方程。
6个方程中任取3个组成方程组，共有C63=20C_6^3=20C63=20种情况，不妨选取以下方程组求解球心：
Mp=N(3)Mp=N \tag 3 Mp=N(3)
其中，
M=[x2−x1y2−y1z2−z1x3−x1y3−y1z4−z1x4−x1y4−y1z4−z1](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=[[(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)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)
经计算，矩阵MMM的行列式值满足：
∣M∣=([x2−x1y2−y1z2−z1]×[x3−x1y3−y1z3−z1])⋅[x4−x1y4−y1z4−z1]=(AB⃗×AC⃗)⋅AD⃗(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,DA,B,C,DA,B,C,D不共面时，混合积(AB⃗×AC⃗)⋅AD⃗(\vec{AB}\times\vec{AC})\cdot\vec{AD}(AB×AC)⋅AD为四面体ABCDABCDABCD的有向体积，不为零，因此矩阵MMM可逆，矩阵方程(3)(3)(3)有唯一解。

三、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;
elsea11 = 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: 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 allpos = 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);
if sum(sqrt((pos(1, :) - sphereCenter(1)).^2 + (pos(2, :) - sphereCenter(2)).^2 + (pos(3, :) - sphereCenter(3)).^2) > radius + 1.0e-6) > 0disp('至少有一个点在球面以外')
endfigure('color', 'w')
draw_sphere(sphereCenter, radius)
hold on
plot3(pos(1, :), pos(2, :), pos(3, :), 'r+')
axis equal tight

任意四面体的外接球/三维空间不共面四点确定唯一球面(附完整代码)相关推荐

任意四面体的外接球的半径(克列尔(A.L.Crelle)公式)
[问题提出]克列尔(A.L.Crelle)公式对任意四面体ABCDABCDABCD,其体积VVV和外接球半径R" role="presentation">RRR满 ...
三角形外接球万能公式_任意四面体的外接球的半径(克列尔(A.L.Crelle)公式)
[问题提出]克列尔(A.L.Crelle)公式对任意四面体$ABCD$,其体积$V$和外接球半径$R$满足$$6RV=\sqrt{p(p-aa_1)(p-bb_1)(p-cc_1)}.$$ 其中$p ...
C++求任意多项式的所有可能的近似根durand kerner roots(附完整源码)
C++求任意多项式的所有可能的近似根durand kerner roots算法 C++求任意多项式的所有可能的近似根durand kerner roots算法完整源码(定义,实现,main函数测试) ...
根据空间中不共面的四个点坐标，求构成任意四面体的内外球
海伦公式: 四面体体积公式六条边分别为a,b,c,a1,b1,c1. a,b,c,a1,b1,c1,其中a与a1,b与b1,c与c1互为对边,那么有三棱锥(四面体)的体积公式为: V=1/12sqr ...
PE结构-空白区手动添加任意代码（附实例代码）
PE之添加任意代码到空白区预备知识 1.查找本机MessageBoxA地址 1.打开OD调试工具拖入要添加的exe程序. 2.在命令中输入 : (输入后按下回车键) 3.点击断点页面即可看到Mess ...
C++：生成给定范围内的所有多维索引。模拟任意数量的嵌套循环的行为(附完整源码)
C++: 生成给定范围内的所有多维索引. 模拟任意数量的嵌套循环的行为 test.hpp test.cpp test.hpp void revers ( int ivec[], int kdim ); ...
输出任意边长的菱形————C语言实践应用（1）（完整源码）
经过一段时间的学习后,想必大家都已经开始摩拳擦掌,迫不及待地想用C语言写一些程序了. 那么今天,我们就来学习C语言中常见的例子--输出任意边长的菱形 ...
最小包围球(附完整代码)
文章目录一.问题描述二.算法步骤三.MATLAB代码一.问题描述给定空间nnn个点,计算最小包围球,使得所有给定点均在球面以内(包括在球面上). 二.算法步骤最小包围球算法可以在最 ...
学点基本功：机器学习常用损失函数小结
(图片付费下载自视觉中国) 作者 | 王桂波转载自知乎用户王桂波 [导读]机器学习中的监督学习本质上是给定一系列训练样本 ,尝试学习的映射关系,使得给定一个 ,即便这个不在训练样本中,也能够得 ...

任意四面体的外接球/三维空间不共面四点确定唯一球面(附完整代码)

文章目录

一、问题描述

二、推导步骤

三、MATLAB代码

任意四面体的外接球/三维空间不共面四点确定唯一球面(附完整代码)相关推荐

最新文章

热门文章