理论参考：
代码实现：
求圆心、半径及方向向量代码
function [cp,R,a,b]=CircleFun(p1,p2,p3)
x1=p1(1);y1=p1(2);z1=p1(3);
x2=p2(1);y2=p2(2);z2=p2(3);
x3=p3(1);y3=p3(2);z3=p3(3);
a1 = (y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2);
b1 = -(x1*z2 - x2*z1 - x1*z3 + x3*z1 + x2*z3 - x3*z2);
c1 = (x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2);
d1 = -(x1*y2*z3 - x1*y3*z2 - x2*y1*z3 + x2*y3*z1 + x3*y1*z2 - x3*y2*z1);
a2 = 2 * (x2 - x1);
b2 = 2 * (y2 - y1);
c2 = 2 * (z2 - z1);
d2 = x1*x1 + y1*y1 + z1*z1 - x2*x2 - y2*y2 - z2*z2;
a3 = 2 * (x3 - x1);
b3 = 2 * (y3 - y1);
c3 = 2 * (z3 - z1);
d3 = x1*x1 + y1*y1 + z1*z1 - x3*x3 - y3*y3 - z3*z3;
M=[a1 b1 c1 0;
a2 b2 c2 0;
a3 b3 c3 0;
0 0 0 1];
cp=inv(M)*[-d1;-d2;-d3;1];
R=sqrt((x1-cp(1))^2+(y1-cp(2))^2+(z1-cp(3))^2);
n=[a1,b1,c1];
n=n/norm(n);%z-vec
a=[x1-cp(1),y1-cp(2),z1-cp(3)];
a=a/norm(a);%x-vec
b=cross(n,a);%y-vec
end参数方程代码
function xyz=generateCircle(cp,R,a,b,num)
xyz=[];
i=linspace(0,2*pi,num);
x=cp(1) + R*cos(i)*a(1) +R*sin(i)*b(1);
y=cp(2) + R*cos(i)*a(2) +R*sin(i)*b(2);
z=cp(3) + R*cos(i)*a(3) +R*sin(i)*b(3);
xyz=[x',y',z'];
end测试代码
clc;
clear all;
close all;
p1=[1,0.5,2]+2*rand(1,3);
p2=[3,-2,4]+3*rand(1,3);
p3=[0,-1.1,2.2]+1.6*rand(1,3);
[cp,R,a,b]=CircleFun(p1,p2,p3);
% [cp1,R1,a1,b1]=CircleFun1(p1,p2,p3);
xyz=generateCircle(cp,R,a,b,100);
figure();
p=[p1;p2;p3];
plot3(p(:,1),p(:,2),p(:,3),'ro');hold on;
plot3(xyz(:,1),xyz(:,2),xyz(:,3));
