用matlab咋三维坐标系内拟合椭球公式

用matlab咋三维坐标系内拟合椭球公式
知道三维坐标系内的一系列的点的坐标,也知道这些点的分布是一个椭球形,怎么用matlab把这个椭球形公式拟合出来?最好是有一段编号的程序,

function my_fit()
% 二维非线性拟合
% 直接将该代码复制到 m文件运行就可以了
% 请仔细看注释,注释写的很清楚
% step0：生成用于拟合的数据

%（以椭球为例,仅为测试,如果有现成数据,请替换此步中 x,y,z 值）
a = 3; %% 方程：x^2/a^2 + y^2/b^2 + z^2/c^2 = 1
b = 4; %% 从而,z = c*sqrt(1 - x^2/a^2 - y^2/b^2)
c = 5; %% 用上半球数据作为待拟合数据

x = -a:0.1:a; %% x,y取值范围
y = -b:0.1:b;
[X, Y] = meshgrid(x,y); %% 生成一个二维的取值范围
[M, N] = size(X);

x = reshape(X, M*N, 1); %% 把矩阵转化为向量
y = reshape(Y, M*N, 1);

p = ((1 - x.^2/a^2 - y.^2/b^2) >= 0); %% 将大于等于0的数值取出（只有这部分才有意义）

x = x(p); %% 生成的值均在上椭球面,如果有现成数据,请将 step0去掉
y = y(p); %% 并直接给 x,y,z赋值
z = c*sqrt(1 - x.^2/a^2 - y.^2/b^2);

% step1：开始拟合,k表示拟合系数,行向量

% 待拟合方程：F = z^2 = c^2 - c^2*x^2/a^2 - c^2*y^2/b^2
% x,y,z 均要先转化为列向量!
% 先把 z 值平方,再进行拟合,很重要!
% 令 c^2 = k(1),c^2/a^2 = k(2), c^2*y^2/b^2 = k(3)
% 求出 k 即得到椭球方程

xdata = [x,y]; %% 将 x,y 数据按列组合到 xdata
ydata = z.^2; %% 先把 z 值平方,再进行拟合
k0 = [1 1 1]; %% k 的运行初值,不会影响最终结果

F = @(k,xdata)k(1) - k(2)*xdata(:,1).^2 -k(3)*xdata(:,2).^2; %% 这句话是拟合函数
[k,resnorm]=lsqcurvefit(F,k0,xdata,ydata); %% 这句话是拟合关键!


% step2：椭圆参数求解
% 根据c^2 = k(1),c^2/a^2 = k(2), c^2*y^2/b^2 = k(3)

c = sqrt(k(1));
a = c/sqrt(k(2));
b = c/sqrt(k(3));

disp('a轴:');
disp(a);
disp('b轴:');
disp(b);
disp('c轴:');
disp(c);

end
再问： 首先谢谢您的帮助，可是椭球的球心不在坐标原点，应该是（x-x1）^2/a^2 + （y-y1）^2/b^2 + （z-z1）^2/c^2 = 1的形式，应该怎么做？而且知道的点的坐标是一系列的不规则的离散点。 不是这种形式x = -a:0.1:a; y = -b:0.1:b。
再答： 1 不规则离散点不会影响结果 2 球心不在原点，方法如下： (x-x1)^2/a^2 +(y-y1)^2/b^2 +(z-z1)^2/c^2 =1 化为： x^2/a^2-2*x1*x/a^2+x1^2/a^2+y^2/b^2-2*y1*y/b^2+y1^2/b^2+z^2/c^2-2*z1*z/c^2+z1^2/c^2 = 1 把z^2项放到一边，代码：（有一行太长这个放不下，注意修改） function my_fit_new() % 日期：2011年12月28日 % 作者：半人马alpha % 由于字数限制，原来有的注释不再写 % 适用于你说的情况 % step0：生成拟合数据（例） a = 3; b = 4; c = 5; x = -a:0.1:a; y = -b:0.1:b; [X, Y] = meshgrid(x,y); [M, N] = size(X); x = reshape(X, M*N, 1); y = reshape(Y, M*N, 1); p = ((1 - x.^2/a^2 - y.^2/b^2) >= 0); x = x(p); y = y(p); z = c*sqrt(1 - x.^2/a^2 - y.^2/b^2); % step1：拟合，k表示系数，行向量 % 待拟合方程：F = z^2 = (-c^2/a^2*x^2) + (c^2/a^2*2*x1*x) + (- c^2/b^2*y^2) + % (c^2/b^2*2*y1*y) + (2*z1*z) + % (-c^2/a^2*x1^2 - c^2/b^2*y1^2 - z1^2 + c^2) % x,y,z 均要先转化为列向量 % k(1) = -c^2/a^2 由k值就可求出椭圆所有参数！！！ % k(2) = c^2/a^2*2*x1 % k(3) = - c^2/b^2 % k(4) = c^2/b^2*2*y1 % k(5) = 2*z1 % k(6) = -c^2/a^2*x1^2 - c^2/b^2*y1^2 - z1^2 + c^2 xdata = [x,y,z]; ydata = z.^2; %% 先把 z 值平方，再进行拟合 k0 = ones(1,6); %% k 的运行初值，不会影响最终结果 F = @(k,xdata) k(1)*xdata(:,1).^2 + k(2)*xdata(:,1) + k(3)*xdata(:,2).^2 + k(4)*xdata(:,2) + k(5)*xdata(:,3) + k(6); [k,resnorm]=lsqcurvefit(F,k0,xdata,ydata); % step2：椭圆参数求解 x1 = -k(2)/k(1)/2; y1 = -k(4)/k(3)/2; z1 = k(5)/2; c = sqrt((z1^2 + k(6))/(1 - 1/a^2*x1^2 - 1/b^2*y1^2)); a = c/sqrt(-k(1)); b = c/sqrt(-k(3)); disp('x1:'); disp(x1); disp('y1:'); disp(y1); disp('z1:'); disp(z1); disp('a轴:'); disp(a); disp('b轴:'); disp(b); disp('c轴:'); disp(c); end
再问： 需要拟合的点的坐标为（0，-174.802,990.048），（0.472，-171.284,995.463），（0.413，-168.639,1003.55），（0.064，-167.862,1019.55），（0，-170.357,1035.44），（0，-172.142,1044.78），（0.215，-174.759,1047.84），（0.171，-176.586,1048.13），（0，-179.832,1043.34），（0,181.589,1040.11），（0，-182.76,1032.62），（0，-184.13,1017.55），（0.113，-183.445,1003.17），用些点拟合成半个椭球。谢谢
再答： function my_fit_new() % 日期：2011年12月29日 % 作者：半人马alpha % 适用于你说的情况 % 你的数据拟合结果是一个旋转双曲面(a,c均为虚数，即 a^2