阿雪的学习记录|MATLAB拟合三维离散点最优空间圆

笑死，宇宙的尽头是matlab和Microsoft Office
本文程序来自这篇博客，以防万一……把我的程序贴上来吧

clear all;
clc;
close all;
%%%   含误差空间圆拟合点  %%%
M=importdata('E:\sjj\0601o\o.txt'); 		%这是我的离散点数据，n行3列
[num dim]=size(M);
 
L1=ones(num,1);
A=inv(M'*M)*M'*L1;       % 求解平面法向量
 
B=zeros(num-1,3);
for i=1:num-1
    B(i,:)=M(i+1,:)-M(i,:);
end
 
L2=zeros(num-1,1);
for i=1:num-1
    L2(i)=(M(i+1,1)^2+M(i+1,2)^2+M(i+1,3)^2-M(i,1)^2-M(i,2)^2-M(i,3)^2)/2;
end
 
D=[B'*B;A']           % 将平面作为限制条件加入方程组
L3=[B'*L2;1]
C=inv(D'*D)*D'*L3   % 求解空间圆圆心坐标
radius=0;
for i=1:num
    tmp=M(i,:)-C';
    radius=radius+sqrt(tmp(1)^2+tmp(2)^2+tmp(3)^2);
end
r=radius/num            %  空间圆拟合半径
 
h1=plot3(M(:,1),M(:,2),M(:,3),'*');
% set(gca,'xlim',[11.4 11.7]);
 
 
%%%%   绘制空间圆  %%%%
n=A;
c=C;
 
theta=(0:2*pi/100:2*pi)';    %  theta角从0到2*pi
a=cross(n,[1 0 0]);          %  n与i叉乘，求取a向量
if ~any(a)                   %  如果a为零向量，将n与j叉乘
    a=cross(n,[0 1 0]);
end
b=cross(n,a);      % 求取b向量
a=a/norm(a);       % 单位化a向量
b=b/norm(b);       % 单位化b向量
c1=(c(1)-10)*ones(size(theta,1),1);
c2=(c(2)-15)*ones(size(theta,1),1);
c3=c(3)*ones(size(theta,1),1);
x=c1+r*a(1)*cos(theta)+r*b(1)*sin(theta);  % 圆上各点的x坐标
y=c2+r*a(2)*cos(theta)+r*b(2)*sin(theta);  % 圆上各点的y坐标
z=c3+r*a(3)*cos(theta)+r*b(3)*sin(theta);  % 圆上各点的z坐标
hold on;
h2=plot3(x,y,z,'-r');
xlabel('x轴')
ylabel('y轴')
zlabel('z轴')
legend([h1 h2],'控制点','拟合圆');
grid on

效果图还不戳，但是要放到论文里的图片能不能发上来啊……就还是不发了，反正效果图还不戳。
然后其实我还需要算拟合率什么的，还妹有找到拟合率的算法，对空间点进行操作真的输麻了，所有事情只有在二维的情况下好算，一涉及到三维，就……毕竟比如说对一堆二维离散点算它的拟合误差，就只需要遍历x坐标，带入拟合曲线的多项式就好了啊，因为是二元方程所以y只会返回一个值，可我现在要把x和y坐标带入圆方程算z，那这合适嘛，这块好像有三种排列组合：(x,y)→z，(x,z)→y，(y,z)→x，而且想想也能感觉到三种方法算误差会有很大差别，而且都不能代表它的实际误差

哦豁要不然算算每组坐标到拟合圆的圆心的距离呢！好像星！
灵感还真是靠胡思乱想才想得出啊……一些野法子，我觉得我要是能早生几百年，我也能是伟大的科学家数学家啥的，然而现在就是我想到一个我觉着好用的方法，一查发现早就有人提出了之类的……但反正殊途同归，对我而言比起灵感乍现，看文献找方法真的太慢了……

每天学一点点，就能少菜一点点，加油( •̀ ω •́ )y