MATLAB 四点定球及三点定圆(完整代码)

哈哈哈看到一个很有意思的想法，顺手写了一下：
怎样快速求圆形或球形的隐函数方程、中心、半径
先从四点定球开始：

四点定球

实际上是看到李扬老师的视频哈哈哈

求通过不共面四点M1(x1,y1,z1)M_{1}\left(x_{1}, y_{1},z_{1}\right)M1(x1,y1,z1),M2(x2,y2,z2)M_{2}\left(x_{2}, y_{2}, z_{2}\right)M2(x2,y2,z2), M3(x3,y3,z3)M_{3}\left(x_{3}, y_{3}, z_{3}\right)M3(x3,y3,z3), M4(x4,y4,z4)M_{4}\left(x_{4}, y_{4}, z_{4}\right)M4(x4,y4,z4)的球面方程

我们设球面方程为：
a(x2+y2+z2)+bx+cy+dz+e=0a\left(x^{2}+y^{2}+z^{2}\right)+b x+c y+d z+e=0 a(x2+y2+z2)+bx+cy+dz+e=0

那么球面上的任意点及已知四点满足如下方程组：

{a(x2+y2+z2)+bx+cy+dz+e=0a(x12+y12+z12)+bx1+cy1+dz1+e=0a(x22+y22+z22)+bx2+cy2+dz2+e=0a(x32+y32+z32)+bx3+cy3+dz3+e=0a(x42+y42+z42)+bx4+cy4+dz4+e=0\left\{\begin{array}{l} a\left(x^{2}+y^{2}+z^{2}\right)+b x+c y+d z+e=0 \\ a\left(x_{1}^{2}+y_{1}^{2}+z_{1}^{2}\right)+b x_{1}+c y_{1}+d z_{1}+e=0 \\ a\left(x_{2}^{2}+y_{2}^{2}+z_{2}^{2}\right)+b x_{2}+c y_{2}+d z_{2}+e=0 \\ a\left(x_{3}^{2}+y_{3}^{2}+z_{3}^{2}\right)+b x_{3}+c y_{3}+d z_{3}+e=0 \\ a\left(x_{4}^{2}+y_{4}^{2}+z_{4}^{2}\right)+b x_{4}+c y_{4}+d z_{4}+e=0 \end{array}\right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧a(x2+y2+z2)+bx+cy+dz+e=0a(x12+y12+z12)+bx1+cy1+dz1+e=0a(x22+y22+z22)+bx2+cy2+dz2+e=0a(x32+y32+z32)+bx3+cy3+dz3+e=0a(x42+y42+z42)+bx4+cy4+dz4+e=0

要使方程组有非零解则要求系数行列式为0，即：
D(x,y,z)=∣x2+y2+z2xyz1x12+y12+z12x1y1z11x22+y22+z22x2y2z21x32+y32+z32x3y3z31x42+y42+z42x4y4z41∣=0D(x, y, z)=\left|\begin{array}{lllll} x^{2}+y^{2}+z^{2} & x & y & z & 1 \\ x_{1}^{2}+y_{1}^{2}+z_{1}^{2} & x_{1} & y_{1} & z_{1} & 1 \\ x_{2}^{2}+y_{2}^{2}+z_{2}^{2} & x_{2} & y_{2} & z_{2} & 1 \\ x_{3}^{2}+y_{3}^{2}+z_{3}^{2} & x_{3} & y_{3} & z_{3} & 1 \\ x_{4}^{2}+y_{4}^{2}+z_{4}^{2} & x_{4} & y_{4} & z_{4} & 1 \end{array}\right|=0 D(x,y,z)=∣∣∣∣∣∣∣∣∣∣x2+y2+z2x12+y12+z12x22+y22+z22x32+y32+z32x42+y42+z42xx1x2x3x4yy1y2y3y4zz1z2z3z411111∣∣∣∣∣∣∣∣∣∣=0

这时候我们发现约束条件D(x,y,z)=0D(x, y, z)=0D(x,y,z)=0是关于x,y,z的隐函数，且由x,y,z的任意性，得出D(x,y,z)=0D(x, y, z)=0D(x,y,z)=0正是我们所要求的球面方程。

当然我们要是想单独求某个系数也非常简单，例如求b只需要求如下行列式:
∣01000x12+y12+z12x1y1z11x22+y22+z22x2y2z21x32+y32+z32x3y3z31x42+y42+z42x4y4z41∣\left|\begin{array}{lllll} 0 & 1 & 0 & 0 & 0 \\ x_{1}^{2}+y_{1}^{2}+z_{1}^{2} & x_{1} & y_{1} & z_{1} & 1 \\ x_{2}^{2}+y_{2}^{2}+z_{2}^{2} & x_{2} & y_{2} & z_{2} & 1 \\ x_{3}^{2}+y_{3}^{2}+z_{3}^{2} & x_{3} & y_{3} & z_{3} & 1 \\ x_{4}^{2}+y_{4}^{2}+z_{4}^{2} & x_{4} & y_{4} & z_{4} & 1 \end{array}\right| ∣∣∣∣∣∣∣∣∣∣0x12+y12+z12x22+y22+z22x32+y32+z32x42+y42+z421x1x2x3x40y1y2y3y40z1z2z3z401111∣∣∣∣∣∣∣∣∣∣

结论成立实践开始，我们编写如下方程：

function [Func,Mu,R]=getBall(X,Y,Z)
syms x y z
symMat=[x.^2+y.^2+z.^2,x,y,z,1];
varMat=[X(1).^2+Y(1).^2+Z(1).^2,X(1),Y(1),Z(1),1;X(2).^2+Y(2).^2+Z(2).^2,X(2),Y(2),Z(2),1;X(3).^2+Y(3).^2+Z(3).^2,X(3),Y(3),Z(3),1;X(4).^2+Y(4).^2+Z(4).^2,X(4),Y(4),Z(4),1];% 计算球隐函数
Func=matlabFunction(det([symMat;varMat]));% 计算各个参数
a=det([1 0 0 0 0;varMat]);
b=det([0 1 0 0 0;varMat]);
c=det([0 0 1 0 0;varMat]);
d=det([0 0 0 1 0;varMat]);
e=det([0 0 0 0 1;varMat]);% 计算球心，半径等信息
Mu=-[b,c,d]./a./2;
R=sqrt(sum(Mu.^2)-e./a);
end

其中返回值Func为球面隐函数，Mu为球心，R为半径

使用实例：

pnt=[0,0,1;1,1,6;1,3,4;9,6,11];
scatter3(pnt(:,1),pnt(:,2),pnt(:,3),'filled')
hold on[~,Mu,R]=getBall(pnt(:,1),pnt(:,2),pnt(:,3));
[X,Y,Z]=sphere(40);
mesh(X.*R+Mu(1),Y.*R+Mu(2),Z.*R+Mu(3),'FaceColor','none')

三点定圆

原理一模一样只不过少了个z而已：

function [Func,Mu,R]=getCircle(X,Y)
syms x y
symMat=[x.^2+y.^2,x,y,1];
varMat=[X(1).^2+Y(1).^2,X(1),Y(1),1;X(2).^2+Y(2).^2,X(2),Y(2),1;X(3).^2+Y(3).^2,X(3),Y(3),1];% 计算圆隐函数
Func=matlabFunction(det([symMat;varMat]));% 计算各个参数
a=det([1 0 0 0;varMat]);
b=det([0 1 0 0;varMat]);
c=det([0 0 1 0;varMat]);
d=det([0 0 0 1;varMat]);% 计算圆心，半径等信息
Mu=-[b,c]./a./2;
R=sqrt(sum(Mu.^2)-d./a);
end

使用实例：

pnt=[1 2;3 4;7 5];
scatter(pnt(:,1),pnt(:,2),'filled')
hold on[~,Mu,R]=getCircle(pnt(:,1),pnt(:,2));
t=linspace(0,2*pi,50);
plot(cos(t).*R+Mu(1),sin(t).*R+Mu(2))

当然这种方法也同样适合很多其他情况，例如两点共线，或者求一些其他曲面的参数之类的，可以自行探索一下。