空间三点画圆代码

时间: 2023-10-04 admin IT培训

空间三点画圆代码

空间三点画圆代码

理论参考:

代码实现:

求圆心、半径及方向向量代码

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-vecend

参数方程代码

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

技术交流邮箱:3531225003@qq