matlab練習程序（橢球擬合）

這次我們來擬合一個橢球，之前也擬合過空間的橢圓，不過當時只用了五個點，方程組應該是欠定的，看看就好。

要擬合橢球，首先設定橢球一般方程：

根據這個方程和已有的空間橢球點數據，利用最小二乘就能得到上面九個參數。

不過有時候我們想要的不是這樣的一般方程，而是橢球的球心和三個半長軸。

下面就來說明如何根據橢球一般方程求取球心和半長軸。

首先把上述方程寫成矩陣形式：

其中xc,yc,zc為球心。

對上式進行展開，得到：

上式和第一個式子相比，求對應一次項相等時的結果。

可以用Mathematica來求：

A = {x - xc, y - yc, z - zc};
B = {{a, 1/2*d, 1/2*e}, {1/2*d, b, 1/2*f}, {1/2*e, 1/2*f, c}};
Factor[Dot[A, B, Transpose[{A}]]]
{a x^2 - 2 a x xc + a xc^2 + d x y - d xc y + b y^2 - d x yc + 
  d xc yc - 2 b y yc + b yc^2 + e x z - e xc z + f y z - f yc z + 
  c z^2 - e x zc + e xc zc - f y zc + f yc zc - 2 c z zc + c zc^2}
CoefficientList[Factor[Dot[A, B, Transpose[{A}]]], {x, y, z}]
{{{{a xc^2 + d xc yc + b yc^2 + e xc zc + f yc zc + c zc^2, -e xc - 
     f yc - 2 c zc, c}, {-d xc - 2 b yc - f zc, f, 0}, {b, 0, 
    0}}, {{-2 a xc - d yc - e zc, e, 0}, {d, 0, 0}, {0, 0, 0}}, {{a, 
    0, 0}, {0, 0, 0}, {0, 0, 0}}}}

得到下面這樣的方程組：

即可求出Xc,Yc,Zc。

半長軸的計算可以參考這里：https://www.zhihu.com/question/47033644/answer/112864757

橢球擬合代碼如下：

clear all;
close all;
clc;

xc = 1.21;
yc = 2.32;
zc = 4.32;

xr = 2.78;
yr = 5.76;
zr = 1.51;

[x,y,z] = ellipsoid(xc,yc,zc,xr,yr,zr,20);

plot3(x(:),y(:),z(:),'.');

x=x(:);
y=y(:);
z=z(:);

X = [x.*x y.*y z.*z x.*y x.*z y.*z x y z];
Y = ones(length(x),1);

C = inv(X'*X)*X'*Y;

M=[C(1) C(4)/2 C(5)/2;
    C(4)/2 C(2) C(6)/2;
    C(5)/2 C(6)/2 C(3)];

Cent = -0.5*[C(7) C(8) C(9)]*inv(M)

S = Cent*M*Cent'+1;
[U,V]=eig(M);

[~,indx] = max(abs(U(1,:)));
[~,indy] = max(abs(U(2,:)));
[~,indz] = max(abs(U(3,:)));

R = [sqrt(S/V(indx,indx)) sqrt(S/V(indy,indy)) sqrt(S/V(indz,indz))]

結果和預設的值是一致的：

參考：https://blog.csdn.net/shenshikexmu/article/details/70143455