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+z12x22+y22+z22x32+y32+z32x42+y42+z42xx1x2x3x4yy1y2y3y4zz1z2z3z411111=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+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])==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')

MATLAB 四点定球及三点定圆(完整代码)_第1张图片

三点定圆

原理一模一样只不过少了个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))

MATLAB 四点定球及三点定圆(完整代码)_第2张图片


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

你可能感兴趣的