过已知4点做球面

碰到一个问题，就是做球面过已知4点，相当于是用这四点确定球的圆心和半径，本想瞄一下算法的，结果还是付费的，干脆自己搞一下。这个问题很简单，球心和半径满足如下关系：

显然可以得到四个方程，分别记为（1），（2），（3），（4），然后用后一个方程减去前面一个方程，可以得到：

其中：

而v就是待求的球心坐标，一旦计算出球心的坐标，则代入到（1）~（4）中任意一个即可求出球的半径。

编写程序如下：

% Q 四行三列，每行表示一个点的x，y，z坐标

function [center,radius]=sphere_4p(Q)
A=zeros(3);
b=zeros(3,1);
for k=1:3
    A(k,:)=2*(Q(k+1,:)-Q(k,:));
    b(k)=Q(k+1,1)^2-Q(k,1)^2+Q(k+1,2)^2-Q(k,2)^2+Q(k+1,3)^2-Q(k,3)^2;
end
if det(A)==0
    error("4 points in a same plane");
end
center=A\b;
radius=0;
for k=1:3
    radius=radius+(Q(4,k)-center(k))^2;
end
radius=sqrt(radius);

%以下为绘图
[x,y,z]=sphere(20);
X=x*radius+center(1);
Y=y*radius+center(2);
Z=z*radius+center(3);
mesh(X,Y,Z);
axis equal
hold on
scatter3(Q(:,1),Q(:,2),Q(:,3),'*');
% 设置坐标轴标签
xlabel('X');
ylabel('Y');
zlabel('Z');

运行该函数：

>> Q=[0 0 0; 3 0 0;0 4 0;0 0 12];
>> [c,r]=sphere_4p(Q)

c =

    1.5000
    2.0000
    6.0000


r =

    6.5000

上面的例子是过（0,0,0），（3,0,0），（0,4,0）和（0,0,12）四点确定球体的中心和坐标，由于取值的特殊性（符合勾股定理），很容易验证其圆心是（3/2,2,13/2）。

得到效果如下图：