用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