#|MATLAB 四点定球及三点定圆(完整代码) 笔记|matlab|线性代数|四点共球

【#|MATLAB 四点定球及三点定圆(完整代码)】哈哈哈看到一个很有意思的想法，顺手写了一下：
怎样快速求圆形或球形的隐函数方程、中心、半径
先从四点定球开始：
四点定球实际上是看到李扬老师的视频哈哈哈
求通过不共面四点 M 1 ( x 1 , y 1 , z 1 ) M_{1}\left(x_{1}, y_{1},z_{1}\right) M1?(x1?,y1?,z1?), M 2 ( x 2 , y 2 , z 2 ) M_{2}\left(x_{2}, y_{2}, z_{2}\right) M2?(x2?,y2?,z2?),M 3 ( x 3 , y 3 , z 3 ) M_{3}\left(x_{3}, y_{3}, z_{3}\right) M3?(x3?,y3?,z3?),M 4 ( x 4 , y 4 , z 4 ) M_{4}\left(x_{4}, y_{4}, z_{4}\right) M4?(x4?,y4?,z4?)的球面方程
我们设球面方程为：
a ( x 2 + y 2 + z 2 ) + b x + c y + d z + e = 0 a\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 ( x 2 + y 2 + z 2 ) + b x + c y + d z + e = 0 a ( x 1 2 + y 1 2 + z 1 2 ) + b x 1 + c y 1 + d z 1 + e = 0 a ( x 2 2 + y 2 2 + z 2 2 ) + b x 2 + c y 2 + d z 2 + e = 0 a ( x 3 2 + y 3 2 + z 3 2 ) + b x 3 + c y 3 + d z 3 + e = 0 a ( x 4 2 + y 4 2 + z 4 2 ) + b x 4 + c y 4 + d z 4 + 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 ) = ∣ 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 ∣ = 0 D(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?+z12?x22?+y22?+z22?x32?+y32?+z32?x42?+y42?+z42??xx1?x2?x3?x4??yy1?y2?y3?y4??zz1?z2?z3?z4??11111?∣∣∣∣∣∣∣∣∣∣?=0
这时候我们发现约束条件 D ( x , y , z ) = 0 D(x, y, z)=0 D(x,y,z)=0是关于x,y,z的隐函数，且由x,y,z的任意性，得出 D ( x , y , z ) = 0 D(x, y, z)=0 D(x,y,z)=0正是我们所要求的球面方程。
当然我们要是想单独求某个系数也非常简单，例如求b只需要求如下行列式:
∣ 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 ∣ \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?+z12?x22?+y22?+z22?x32?+y32?+z32?x42?+y42?+z42??1x1?x2?x3?x4??0y1?y2?y3?y4??0z1?z2?z3?z4??01111?∣∣∣∣∣∣∣∣∣∣?
结论成立实践开始，我们编写如下方程：

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])==0); % 计算各个参数 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 z 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])==0); % 计算各个参数 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))