Space three point circle code

Theoretical reference ：

Code implementation ：

Find the center of a circle 、 Radius and direction vector code

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

Parametric equation code

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

Test code

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));

Technical exchange mailbox ：[email protected]