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

Question

各位大侠，知道三维坐标系内的一系列的点的坐标，也知道这些点的分布是一个椭球形，怎么用matlab把这个椭球形公式拟合出来？最好是有一段编号的程序，谢谢！

Answer 1

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<0,c^2<0)
% 我按拟合出的参数给你把图画了一下，是旋转双曲面的一支

% step0：生成拟合数据（例）
x = [0,0,0,0,0,0,0,0.064,0.113,0.171,0.215,0.413,0.472]';
y = [-174.802,-170.357,-172.142,-179.832,181.589,-182.760,-184.130,-167.862,-183.445,-176.586,-174.759,-168.639,-171.284]';
z = [990.048,1035.44,1044.78,1043.34,1040.11,1032.62,1017.55,1019.55,1003.17,1048.13,1047.84,1003.55,995.463]';

% 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) - k(1)*x1^2 - k(3)*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

追问

大侠，十分感谢，可是还有些技术上的细节问题，能加你qq么？

追答

我的QQ:944096506