Fit ellipsoid to (x,y,z) data

15 次查看（过去 30 天）
Geetartha Dutta
Geetartha Dutta 2023-10-25
评论： Geetartha Dutta 2023-11-1
I have a 3D dataset having (x,y,z) coordinates. The x and y values are equally spaced (regular grid). How can I fit an ellipsoid of the form (x-p)^2/a^2 + (y-q)^2/b^2 + (z-r)^2/c^2 , where (p,q,r) are the coordinates of the center of the ellipsoid, and a,b,c are the radii?
  7 个评论
显示 5更早的评论
Matt J
Matt J 2023-10-26
编辑：Matt J 2023-10-26
I know that there seems to be two modes in the data
Looks like a lot more than that. I can't tell which is supposed to be the "greater" mode. In any case, if you want a good fit in a particular region, you will have to prune the data to exclude the other regions.
Geetartha Dutta
Geetartha Dutta 2023-10-26
Attached is the pruned data. It would be great if I could get a reasonably good fit to this data.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2023-10-26
编辑：Matt J 2023-10-26
Ran in:
I'm finding that a decent fitting strategy is to first fit with a Gaussian, but then use the parameters of the Gaussian to construct an ellipsoid hemisphere. For the Gaussian fitting, I used gaussfitn, which is downloadable from,
https://www.mathworks.com/matlabcentral/fileexchange/69116-gaussfitn
load xyz
[maxval,i]=max(z(:));
mu0=[x(i);y(i)];
D0=min(z(:));
opts={'FunctionTolerance',1e-14, 'OptimalityTolerance',1e-14, 'StepTolerance',1e-14};
G0={D0,maxval-D0,mu0,100*eye(2)};
LB={0,0,[],[]};
UB={D0,maxval,[],[]};
G = gaussfitn([x(:),y(:)],z(:),G0,LB,UB,opts{:});
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
%Disaply surfaces
[Zg,Ze]=getSurf(x,y,G);
surf(x,y,z,'FaceAlpha',0.5,'FaceColor','b');
surface(x,y,Ze,'FaceColor','r'); xlabel X, ylabel Y
legend('Raw Data','Fit')
function [Zg,Ze]=getSurf(x,y,G)
[D,A,mu,sig]=deal(G{:});
sz=size(x);
xy=[x(:),y(:)]'-mu;
Zg=D+A*exp(-0.5*sum( (sig\xy).*xy,1)); Zg=reshape(Zg,sz); %Gaussian Fit
Ze=D+A*sqrt(1-sum( (sig\xy).*xy)); Ze=reshape(Ze,sz); %Ellipsoid Fit
end
  6 个评论
显示 4更早的评论
Matt J
Matt J 2023-10-31
编辑：Matt J 2023-10-31
You should set the complex values to NaN. They correspond to (x,y) outside the footprint of the ellipsoid.

请先登录，再进行评论。

更多回答（2 个）

Torsten
Torsten 2023-10-25
Matt J
Matt J 2023-10-26
编辑：Matt J 2023-10-26
Ran in:
Using quadricFit from,
https://uk.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics
%%%%%%%%%%%Fake input data
[X,Y,Z] = sphere;
[X,Y,Z]=deal(1+40*X, 2+20*Y,3+30*Z); %stretch into an ellipsoid
surf(X,Y,Z); axis equal
%%%%%%%%%%% Do the fit
XYZ=[X(:),Y(:),Z(:)]';
[XYZ,T]=quadricFit.homogNorm(XYZ);
X=XYZ(1,:).';
Y=XYZ(2,:).';
Z=XYZ(3,:).';
e=+ones(size(X,1),1);
M= [X.^2, [], [], X, ...
Y.^2, [], Y, ...
Z.^2 Z, ...
e];
coeffs=quadricFit.mostnull(M);
ABCDEFGHIJ=zeros(1,10);
ABCDEFGHIJ([1,4,5,7:10])=coeffs;
ABCDEFGHIJ=num2cell(ABCDEFGHIJ);
[A,B,C,D,E,F,G,H,I,J]=deal(ABCDEFGHIJ{:});
Q=[A, B, C; %D
0 E, F; %G
0 0 H];%I
%J
Q=Q/2+Q.'/2;
W=T.'*[Q,[D;G;I]/2;[D,G,I]/2,J]*T;
Q=W(1:3,1:3);
x0=-Q\W(1:3,end);
T=eye(4); T(1:3,4)=x0;
W=T.'*W*T; W=-W/W(end);
rad=sqrt(1./diag(W(1:3,1:3)));
[a,b,c]=deal(rad(1),rad(2),rad(3)) %ellipsoid radii
a = 40.0000
b = 20
c = 30.0000
[p,q,r]=deal(x0(1),x0(2),x0(3)) %ellipsoid center coordinates
p = 1.0000
q = 2.0000
r = 3.0000
  2 个评论
Geetartha Dutta
Geetartha Dutta 2023-10-26
I tried the above code using my data, and it gives complex values for a and b. I am not sure why.
Matt J
Matt J 2023-10-26
Attach your xyz data in a .mat file, so it can be examined.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Curve Fitting Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by