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

一、问题描述

  给定三维空间不共面的四个点$A,B,C,D$，求由此四点确定的球面方程$(x-x_0{)}^2+(y-y_0{)}^2+(z-z_0{)}^2=R^2$。

二、推导步骤

方法一

  设坐标$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)$，由于此四点在球面上，因而：
$$\{\begin{array}{ll}(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)\end{array}\text{\hspace{1em}(1)}$$
  式$(1.1)$至$(1.4)$两两做差，整理可得：
$$\begin{array}{rl}& (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\end{array}\text{\hspace{1em}(2)}$$
  其中，$i,j=1,2,3,4$且$i\ne j$，即式$(2)$有$C_4^2=6$个方程。
  6个方程中任取3个组成方程组，共有$C_6^3=20$种情况，不妨选取以下方程组求解球心：
$$Mp=N\text{\hspace{1em}(3)}$$
  其中，
$$M=\left[\begin{array}{ccc}{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{array}\right]\text{\hspace{1em}(4)}$$
$$N=\left[\begin{array}{c}[(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{array}\right]\text{\hspace{1em}(5)}$$
  经计算，矩阵$M$的行列式值满足：
$$\begin{array}{rl}|M|& =(\left[\begin{array}{c}{x}_2-x_1\\ y_2-y_1\\ z_2-z_1\end{array}\right]\times \left[\begin{array}{c}{x}_3-x_1\\ y_3-y_1\\ z_3-z_1\end{array}\right])\cdot \left[\begin{array}{c}{x}_4-x_1\\ y_4-y_1\\ z_4-z_1\end{array}\right]\\ & =(\overrightarrow{AB}\times \overrightarrow{AC})\cdot \overrightarrow{AD}\end{array}\text{\hspace{1em}(6)}$$
  当$A,B,C,D$不共面时，混合积$(\overrightarrow{AB}\times \overrightarrow{AC})\cdot \overrightarrow{AD}$为四面体$ABCD$的有向体积，不为零，因此矩阵$M$可逆，矩阵方程$(3)$有唯一解。

方法二

  将球面方程展开成：
$$x^2+y^2+z^2-2x_{0}x-2y_{0}y-2z_{0}z+x_0^2+y_0^2+z_0^2-R^2=0\text{\hspace{1em}(7)}$$
  根据多项式的项和系数特征可以将式(7)写成：
$$x^2+y^2+z^2+ax+by+cz+d=0\text{\hspace{1em}(8)}$$
  将给定4个点$(x_i,y_i,z_i)(i=1,2,3,4)$带入式(8)，可以得到关于待定系数$a,b,c,d$的四元一次线性方程组。求解四元一次线性方程组，可以求得球心坐标：
$$\{\begin{array}{ll}{x}_0=-a/2& \\ y_0=-b/2& \\ z_0=-c/2& \end{array}\text{\hspace{1em}(9)}$$
  球半径：
$$R=\sqrt{(a/2{)}^2+(b/2{)}^2+(c/2{)}^2-d}\text{\hspace{1em}(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